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ABSTRACT 


A  theoretical  and  experimental  investigation  is  presented  for 
the  problem  of  the  transient  growth  of  a  frozen  shell  of  ice  that  develops 
when  superheated  water  flows  through  a  vertical  cylinder,  whose  wall  is 
cooled  convectively  to  a  temperature  below  the  freezing  temperature.  The 
solution  of  the  theoretical  problem  is  developed  from  the  asymptotic 
solution  for  the  temperature  distribution  in  the  ice,  with  prescribed 
convection  as  the  boundary  condition.  This  analysis  was  broken  into  two 
parts;  a  zone  which  is  ice-free  and  a  zone  in  which  the  ice  formation 
takes  place.  The  experimental  results  obtained  agree  well  with  the  corre¬ 
sponding  theoretical  results  in  the  region  where  the  theory  is  applicable. 
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CHAPTER  I 


INTRODUCTION 

Transient  solidification  of  a  superheated  liquid  flowing  through 
a  tube,  whose  wall  temperature  is  lowered  below  that  of  the  liquid  freezing 
temperature,  has  a  complicated  history.  As  the  liquid  undergoes  a  change 
to  the  solid  phase,  thermal  energy  in  the  form  of  the  latent  heat  of  fusion 
is  released  from  the  moving  solid-liquid  interface  where,  in  addition,  heat 
may  be  removed  from  the  flowing  liquid. 

This  type  of  solidification  problem  arises  in  many  familiar  and 
important  circumstances:  for  example;  a  buried  pipe  in  a  region  of  sporadic 
permafrost,  the  casting  of  metals,  and  freezing  of  foods.  In  some  circum¬ 
stances  the  tube  wall  may  be  cooled  convectively ,  e.g.  a  water  line  exposed 
to  a  "cold"  meteorological  environment.  Another  area  where  problems  of  this 
type  arise  is  in  space  technology.  An  example  is  provided  by  the  hydraulic 
lines  of  space  satellites  which  may  be  subjected  to  extreme  environmental 
temperatures;  if  proper  precautions  are  not  taken  solidification  may  occur. 
Yet  another  instance  is  in  propulsion  devices  that  utilize  liquid-liquid 
heat  exchangers  in  which  one  liquid  is  below  the  freezing  temperature  of 
the  other:  if  the  flow  rate  or  the  temperature  of  the  "warm"  liquid  is 
low  enough,  a  frozen  layer  will  form  on  the  exchanger  walls. 

A  survey  of  the  literature  reveals  that  little  information  is 
available  on  the  solidification  of  a  liquid  flowing  through  a  tube  with  a 
convective  boundary  condition  at  the  moving  interface  or  a  convective 
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boundary  condition  at  the  tube  wall.  Brush  [1]*  discussed  the  principles 
governing  the  freezing  of  water  in  mains.  He  presented  no  theoretical 
or  experimental  data  but  suggested  that  the  formation  of  ice  in  water 
mains  is  dependent  upon  the  temperature  of  the  water  and  its  flow  rate 
as  well  as  on  the  wall  temperature. 

London  and  Seban  [2]  describe  a  general  approximate  method  for 
the  solution  of  the  one-coordinate  freezing  problem  with  application  to 
the  ice  formation  within  various  geometrical  boundaries.  Their  analytical 
solutions  were  obtained  from  the  construction  of  a  thermal  circuit,  con¬ 
sidering  the  sensible  heat  negligible  relative  to  the  latent  heat  of 
fusion  of  the  liquid.  They  presented  graphic  and  algebraic  solutions 
for  ice  formation  at  spherical,  cylindrical  and  plane  boundaries  along 
with  an  investigated  degree  of  approximation  of  these  solutions.  It  was 
shown  that  the  analytical  predictions  of  freezing  times  were  comparable  to 
those  obtained  experimentally. 

Kreith  and  Romie  [3]  presented  solutions  for  determining  the 
position  of  the  solid-liquid  interface  by  considering  the  limited  condition 
of  a  constant  fusion  velocity  for  cylinders  and  spheres  with  the  water 
initially  at  the  fusion  temperature. 

Stephenson  [4]  studied  the  problem  of  water  flowing  through  a 
pipe  exposed  to  an  environment  at  temperatures  below  the  freezing  point. 

He  presented  an  empirical  equation,  using  thermal  resistances,  to  relate 
three  parameters  of  the  problem;  minimum  mass  flow  rate  required  for  ice 


^Numbers  in  brackets  denote  references  on  page  68. 
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not  to  form,  temperature  of  the  water  entering  the  exposed  pipe,  and 
resistance  to  heat  transfer  from  the  water  to  the  environment. 

Zerkle  and  Sunderland  [5]  discussed  the  effect  of  liquid 
solidification  at  the  inner  surface  of  a  horizontal  circular  tube  upon 
steady  laminar  flow  heat  transfer:  they  assumed  steady  state  conditions, 
i.e.  neglected  latent  heat  effects.  They  found  that,  as  could  be  ex¬ 
pected,  the  amount  of  growth  of  the  solid-liquid  interface  is  dependent 
on  the  amount  of  superheat  in  the  liquid  and  its  flow  rate.  Their  experi¬ 
ments  were  carried  out  with  water  in  horizontal  tubes  and  it  was  found 
that  free  convection  produced  considerable  deviations  from  the  theoretical 
results  which  ignored  free  convection. 

Beaubouef  and  Chapman  [6]  analyzed  the  freezing  of  a  solid 

phase  from  a  fluid  flowing  past  a  "cold"  surface.  A  finite-difference 
technique  was  employed  to  determine  the  thickness  of  the  solid  phase  de¬ 
posited  by  the  flowing  liquid,  as  a  function  of  time  and  location  on  the 
surface. 

Savino  and  Siegel  [7]  presented  analytical  predictions  and 
experimental  data  for  the  transient  growth  of  a  frozen  layer  which  forms 
when  a  flowing  warm  liquid  is  in  contact  with  a  plane  wall  that  is  cooled 
convectively  from  the  opposite  side.  Three  analytical  procedures  were 
employed  and  compared  for  accuracy  and  convenience  in  application;  one  of 
these  is  an  iterative  technique.  Experimental  values  of  ice  thickness 
and  plate  surface  temperature  obtained  throughout  the  transient  growth 
period  compared  favorably  with  the  theoretically  predicted  values. 
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Gort  [8]  studied  the  problem  of  freezing  inside  a  long  vertical 
cylinder  with  and  without  a  forced  internal  flow.  The  theoretical  solution 
which  he  used  is  an  approximate  solution  for  arbitrary  departures  of  the 
cylinder  wall  temperature  beneath  the  freezing  temperature,  with  the  fluid 
at  its  freezing  temperature.  Experimental  work  on  this  Dirichlet  problem 
was  conducted  for  no-flow  and  flow  cases  with  particular  emphasis  on  the 
influence  of  flow  (and  superheat)  on  the  ice  formation  studied. 

The  recent  work  of  Des  Ruisseaux  and  Zerkle  [9]  is  an  extension 
of  the  work  of  Zerkle  and  Sunderland  [5]  where  a  theoretical  technique 
is  developed  for  predicting  conditions  under  which  a  hydraulic  system 
freezes  shut.  System  pressure  drop  as  a  function  of  Reynolds  number  at 
steady-state  conditions  was  determined  analytically  for  various  uniform 
wall  temperatures.  The  theoretical  results  indicate  a  minimum  pressure 
drop  which  should  be  maintained  across  a  system  to  prevent  it  from  freezing 
shut. 

•  • 

Very  recently,  Ozisik  and  Mulligan  [14]  presented  an  analytical 
investigation  of  the  transient  freezing  of  a  liquid  flowing  in  a  circular 
tube  under  the  assumption  of  a  constant  wall  temperature  which  is  lower 
than  the  liquid  freezing  temperature:  constant  properties,  a  slug-flow 
velocity  profile  and  quasi-steady  state  heat  conduction  In  the  solid  phase 
were  assumed.  The  method  of  solution  involved  the  determination  of  an  ex¬ 
pression  for  the  dimensionless  temperature  of  both  the  solid  and  liquid 
phases  in  terms  of  the  unknown  solid-liquid  interface  radius.  These  ex¬ 
pressions  are  coupled  through  the  heat  balance  at  the  interface  resulting 
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in  a  single  ordinary  differential  equation  for  the  interface  radius.  The 
analysis  presented  the  variation  of  the  local  heat  flux  and  the  solid- 
liquid  interface  profile  during  freezing  as  a  function  of  time  and  position 
along  the  tube.  The  analysis  produced  steady  state  heat  transfer  rates 
and  axial  profiles  which  were  compared  with  the  results  of  Zerkle  and 
Sunderl and. 

The  purpose  of  this  thesis  is  to  study  the  problem  of  transient 
freezing  in  a  long  vertical  cylinder  with  a  forced  internal  flow.  An 
approximate  theoretical  solution  for  the  interface  history  includes  the 
convective  heat  transfer  condition  imposed  at  the  inward  moving  interface 
of  the  frozen  layer.  The  analysis  is  also  concerned  with  a  convective 
boundary  condition  on  the  outside  of  the  cylinder  wall,  whose  temperature 
is  below  the  freezing  point.  Over  an  initial  length  of  the  cylinder  the 
cooling  is  insufficient  to  cause  solidification  and  therefore  an  ice- 
free  zone  exists.  From  this  point  on,  ice  grows  radially  inwards  at  a 
rate  which  is  found  to  increase  with  the  distance  along  the  tube.  The 
theoretical  results  are  used  to  show  how  the  ice  thickness  and  pressure 
drop  vary  with  length. 

Representative  experimental  results  are  compared  with  the 
theoretical  predictions  and  indicate  the  accuracy  and  range  of  validity 
of  the  theory. 
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PART  I 


THEORETICAL  ANALYSIS 


CHAPTER  II 


SOLUTION  OF  THE  ICE-FREE  PROBLEM 


2.1  Governing  Equations 

Consider  superheated  water  flowing  vertically  upwards  in  steady 
laminar  motion,  occupying  the  cylindrical  domain  shown  in  Figure  2.1. 

The  water  is  assumed  to  have  a  uniform  temperature  and  a  fully  developed 
velocity  profile  at  the  thermal  entrance,  i.e.  z  =  0.  A  prescribed  con¬ 
vective  boundary  condition  is  applied  at  the  wall  of  the  cylinder,  above 
z  =  0,  for  times  greater  than  zero.  The  tube  wall  has  negligible  thermal 
resistance  in  the  cooling  region. 

The  approximate  form  of  the  differential  equation  governing  heat 
conduction  in  the  water  downstream  from  the  thermal  entrance,  with  radial 
symmetry  and  no  sources,  is  given  by 


9T 


!!li  +  i!V:L(v  _ 

3r2  r  3r  “*  2  32 


£ 


9T 


£' 


+  vr  37*) 


2.1-1 


Boundary  conditions  (including  convective  cooling)  may  be  written 


9T 


£> 


-  0 

3  r=0 


k*  (5T>  ■  h0(Tw  -  Tc> 
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DIRECTION  OF  FLOW 


FIGURE  2.1  COORDINATE  SYSTEM  OF  A  TUBE  WITH  SOLIDIFICATION 
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T£(r,°)  -  Tq 


(An  order  of  magnitude  analysis  -  see  Appendix  B  -  reveals 
in  the  axial  direction  may  be  neglected  since  p'  ^  «  1). 

With  the  velocity  profile  assumed  fully  developed,  i.e.  v^ 
cable  equation  is 


that  conduction 


=  0,  the  appli- 


32T 


3r 


i  +  1 

r  3r 


vz  3T* 


°* 3z 


Normalizing  equation  2.1-2  gives 


where 


fV  i_!!i 

3r+2  r+  8r+ 


+  3$ 

u 


£ 


r0Ped 


V  = 


S_ 

2 


2.1-2 


2.1-3 


For  hydrodynamical ly  fully  developed  flow  the  parabolic  velocity 


v  =  2V[1  -  (f-)2] 
z  r0 


profile  is  applicable 


01  -  (0.-r)sT 
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Normalizing  yields 


u+  =  2(1  -  r+2) 


2.1-4 


Substituting  into  equation  2.1-2,  the  final  form  of  the 
differential  equation  is 


3r+2  r+  3 r+  "  ”  r  }  ^7 


2.1-5 


The  boundary  conditions  then  are 


(~f)  +  "  0 

9r  r  =0 


8$ 

(-4)  .  =  -  Nu  $ 

'  +'  +  ,  c  w 

9r  r  =1 


$£(r  ,0)  =  1 


where 


Nu_  = 


roho 


c  k 


l 


$  = 


T  -  T 
w  c 


w  T0-Tc 


Equation  2.1-5  and  the  succeeding  boundary  conditions  are  linear 
and  homogeneous,  and  as  such  can  be  solved  using  the  method  of  separation 
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of  variables  by  taking 


Vr+,z+)  =  R(r+)-Z(z+)  •  2.1-6 

This  leads  to  two  ordinary  differential  equations^ 

v  +  \h  =  0  2.1-7 

R"  +  -L  R'  +  x2r(1  _  r+2}  =  0  _  2.1-8 

r 

The  first  of  these  is  satisfied  by  an  exponential  function, 
i.e.  Z  =  C  exp(-  X  z  ).  The  second  equation  is  of  the  Sturm-Liouvi 1 le 
type,  where  the  solutions  that  satisfy  the  boundary  conditions  can  be 
found  for  an  infinite  number  of  eigenvalues  X.  That  is 

R(r+)  =  AR1(r+)  +  BR2(r+). 


Yielding 


$£(r+,z+)  =  [AR-,  (r+)  +  BR2(r+)]  C  exp(-  X2z+)  . 
Application  of  the  boundary  conditions  leads  to 


B  =  0 


t,  . 

3 


first  and  second  derivatives  with  respect  to  the  argument 
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Nu 


1  Vr+> 

+ 


c  R-j  '  (r+)  r+=l 


=  0 


2.1-9 


Thus  the  solution  of  equation  2,1-5  may  be  written  in  the  form 


GO 

$£(r+,z+)  =  l  cn  Rp(r+)  exp  (-  An2z+) 


2.1-10 


n=0 


where  A n  are  the  roots  of  the  transcendental  equation  and  Rn(r+)  are 
the  eigenfunctions  corresponding  to  the  function  R-|(r+).  Applying  the 
initial  condition  the  coefficients  cn  are  determined  by  the  relation  [15] 

1  r  (1-r  )  R(r  )  dr 

0  n 

cn  (-1  +/-,  +2\  D  2/  +\  ,  +  ^ 

/  r  (1-r  )  Rn  (r  )  dr 

0  n 

since  R  (r+)  is  an  orthonormal  function, 
n 

Following  the  work  of  Sellars,  Tribus  and  Klein  [10],  the  solu¬ 
tion  of  equation  2,1-8  is  sought  for  higher  eigenvalues.  This  was  ac¬ 
complished  by  obtaining  an  asymptotic  solution  for  equation  2.1-8  for 
An  00  .  The  resulting  equation  for  R(r+)  is  the  WKB  approximation  and 
is  valid  for  0  <  r+  <  1  and  large  A.  Examination  of  equation  2.1-8  shows 
that  when  r  is  very  small,  i.e.  A  (1  -  r  )  •+  A  ,  the  classical  solution 
behaves  as  a  Bessel  function  of  zero  order  of  the  first  kind:  thus  R(r+)  = 
Jg(Ar  ).  For  large  Ar  ,  the  asymptotic  solution  for  JQ ( Ar  )  is  matched  to 
the  WKB  approximation  to  determine  the  coefficients  of  the  latter  for  r+ 


ffiJnsbneoanfi'i}  erli  aiocn  9rid  s'ic  X  9^riw 
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small.  Since  the  resulting  solution  has  a  singularity  at  r+  1 ,  an 
alternate  solution  is  required  because  a  boundary  condition  is  to  be 
imposed  at  r+  =  1.  This  is  accomplished  by  making  a  change  of  variable, 
i.e.  x+  =  1  -  r+,  and  substituting  into  equation  2.1-8  for  small  x+. 

Yet  another  substitution  of  a  new  variable  yields,  for  large  A,  a  differ¬ 
ential  equation  which  results  in  a  further  solution  to  equation  2.1-8. 

This  new  solution  is  matched  to  the  previously  obtained  solution  to  evaluate 
the  formers'  unknown  constants.  The  above  procedure  thus  provides  a  solution 
valid  for  small  x+  (near  the  wall).  The  solution  is  [10] 

+3/2 

n / .,+  \  _  2  /757#+  r Att  tt ^  ,  /A/8  x  \ 

R(x  )  ~  ^  / 2x  L s i n (  4  “  3)  ^1  (  2  * 

3 

.  / Att  2i\\  ,  /A/8  x+3/^n  o  i  no 

-  sm(-^ - 3-)  J1  ( - 3 - )J  2.1-12 

"I 


2.2  Eigenvalues 


R(r+) 

From  equation  2.1-12  the  ratio  — * 


is  found  to  be 


(r  )r  =1 


R(r+) 


RV) 


=  R(x+)/  - 
r+=l  8x 


x+=0 


q3  r/4x  •  (X £  2tt_\ 

3  1(3)  sin  -  3) 

1  2 

23  r(|)  X3  sin(^-  -  |) 


2.2-1 


This  may  now  be  substituted  into  the  convective  boundary  condition 
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(equation  2.1-9)  resulting  in  the  transcendental  equation. 


2_ 

33  r(f)  sin  -  §4 


1  2 

23  r(f)  X3  sin  (^  -  f) 


=  0 


2.2-2 


This  equation  leads  to  an  infinite  set  of  eigenvalues,  each  set  depending 
on  the  choice  of  Nusselt  number  Nu  .  As  Nu  -*  °°  ,  the  eigenvalues  tend 
to  those  obtained  by  Sellars,  et  al  [10]  since  the  boundary  condition  then 
becomes  a  Dirichlet  type,  i.e.  a  constant  wall  temperature,  in  which  case 
it  is  evident  that  A  =  4n  +  ^-  . 

J 


2.3  Ice-Free  Length 

Once  the  eigenvalues  have  been  determined  from  the  boundary 
condition  at  the  wall,  the  length  of  the  ice-free  zone,  z*  ,  may  be 
determined.  This  is  found  from  equation  2.1-10  such  that 


V1’  2e+)  =  4e 


At  this  point,  where  the  wall  temperature  is  at  the  freezing  temperature, 
the  corresponding  temperature  profile  must  be  calculated  again  using 
equation  2.1-10.  This  temperature  profile  will  become  the  input  profile 
to  the  freezing  zone. 

The  theoretical  results  are  shown  graphically  in  Figure  2.2 
which  demonstrates  the  variation  in  the  length  of  the  ice-free  zone 
versus  superheat  ratio  [a  =  (Tq  -  T^)/(T^  -  Tc)]  for  various  Biot 
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FIGURE  2.2  VARIATION  OF  z  WITH  SUPERHEAT  RATIO 
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numbers  [Bi  =  (k  /k.)  Nu  ].  As  shown,  the  ice-free  length  decreases 

v  X/  I  L 

with  decreased  superheat  or  increased  Biot  number  for  a  particular 
Reynolds  number.  This  verifies  what  was  expected  a  priori. 


CHAPTER  III 


SOLUTION  OF  THE  FREEZING  PROBLEM 
3.1  Governing  Equations 

This  chapter  is  concerned  with  the  transient  freezing  of  water 
flowing  vertically  upwards  in  steady  laminar  motion  through  a  cylindrical 
tube  with  solidification  at  the  inner  wall  surface.  Consider  the  flowing 
water  and  its  corresponding  ice  shell  occupying  the  cylindrical  domain 
shown  in  Figure  2.1.  The  superheated  water  now  has  a  modified  temperature 
profile  at  the  leading  edge  of  the  freezing  zone  due  to  cooling  in  the 
ice-free  zone.  Since  the  water  is  cooled  further  as  it  flows  through  the 
freezing  zone,  its  mean  temperature  approaches  the  freezing  temperature 
and  the  thickness  of  the  ice  shell  increases.  Due  to  this  variation  in 
the  cross-section  a  two-dimensional  velocity  distribution  is  produced. 

Assuming  the  ice  properties  constant,  the  conduction  equation 
within  the  ice  layer,  with  radial  symmetry  and  no  sources,  is  given  by 


,  3Ti 

ar2 


r  Sr 


a.  St 


3.1-1 


Where  axial  conduction  has  been  shown  by  an  order  of  magnitude  analysis 
(Appendix  B)  to  be  negligible.  The  system  is  subject  to  the  following 
initial  and  boundary  conditions 
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Ti(6,t)  =  T  f 


3T . 

k.(^)  =  hn(T  -  T  ) 

i  v3r  '  Cr  w  c' 

r=r0 


W0)  ■  Tf 


The  temperature  distribution  in  the  liquid  downstream  from 
the  entrance  to  the  freezing  zone  is  described  by  the  energy  equation 


9T£  rl  A.  (  i 

Vr  3r  Vz  3z  a£^r  3r  3r 


3,1-2 


With  the  corresponding  boundary  conditions 


Vr’0)  =  T0 


T£(6,z)  Tf 


3T 


£' 


^  =  ° 
r=0 


Conservation  of  mass  is  represented  by 


3(rv  )  Sv 

1 - —  +  — -  =  o 

r  9r  3z 


3,1-3 


or 


2tt  /  v  r  dr  =  Q 
0  z 
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Taking  as  an  approximation 


vz(r,z) 


2V(^-)2  [1 


<£)  ] 


3.1-4 


and  using  equations  3,1-3 


vr(r,z) 


r  r 


=  2  V 


0  jd6 
dz 


Q 

\  'iV 


[1  - 


(j)  ]• 


3.1-5 


These  equations  imply  that  the  flow  cross  section  is  slowly  varying  and 
vr  d6 

reveal  that  --  <<  1  i f  ~  «  1 ,  After  substitution  of  the  velocity 
z 

components,  equation  3,1-2  becomes 


rn  2  2  ,,  9T0  9T  9T0 


3,1-6 


Following  Zerkle  and  Sunderland  [5],  the  above  equation  is 
normalized,  yielding 


6* 


J_  n  -  f— dd*  — -  +  =  r — —  +  J. — -^\ 


924) 


9r* 


9<J>, 

2  1  r*  9r* 


3.1-7 


<J>A(r*,0)  =  1 


4»a(6*,z*)  =  0 


with  the  boundary  conditions 
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34>, 

( _ h 

[dr*J 


=  0 


r*=0 


where 


<t>4(r*,z*)  = 


0 


l 


6ci  ’  0£  T{'"  Tf  ’  ^  T°  "  Tf 


r*  =  -,{*  =  -,  z*  =  — | — 
r0  r0  r0Ped 


It  should  be  emphasized  at  this  point  that  the  axial  position  is  now 

measured  from  the  leading  edge  of  the  freezing  zone, 

r 

When  the  substitution,  n  =  j  ,  is  used  in  equation  3.1-7,  the 
resulting  equation  and  boundary  conditions  are  of  the  same  form  as  the 
system  describing  the  classical  Graetz  problem 


O  -n2)^  =  [ 


dr\ 


Z  1 

2  n  9n  J 


3.1-8 


and  the  corresponding  boundary  conditions 


4>£(n,0)  =  1 


4>a(1.z*)  =  0 


3<J> 


Z' 


^  n  =  ° 

1  n=0 


The  non-dimensional  temperature  distribution  is  then  given  by 


cj)  (n,z)  =  l  c.  R.(n)  exp  (-  X.  z*) 

X>  j  _  Q  J  J  J 


3.1-9 
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with  the  temperature  gradient  at  the  ice-water  interface  given  by 


[_  ^ 


oo 


3n  0>z*)]=  2  [-  CjRj  '  ( 1  )/2]  exp  (-x/z*) 


3.1-10 


and  the  heat  transfer  by 


q  = 


2MW 


00 


j=0 


c . 
2 


R . 1  (1)  exp  (-A.  z*) 


3.1-11 


In  the  present  problem,  however,  the  values  of  the  coefficient, 
c.,  cannot  be  taken  from  Sellars, et  al  [10]  as  done  by  Zerkle  and 

J 

Sunderland  [5].  This  is  due  to  the  fact  that  the  temperature  profile  at 
the  leading  edge  of  the  freezing  zone  is  not  a  block  profile  but  has 
been  modified  by  convection  in  the  ice-free  zone.  Thus  we  must  write  [15] 


/%  (n,0)n(l-n2)Mn)dn 
0  J 


J^oO-o2)  R,2  (n )  dn 
0  J 


3.1-12 


where 


oo 


<J>A(n.O)  = 


<W  +  (VV  L  cn  exp  (-X„V> 

_  n-u _ _ 

T0  "  Tf 


+  + 
noting  that  at  zg  ,  n  is  identically  r  . 

Consider  now  the  heat  balance  at  the  interface.  The  heat  flux 
through  the  frozen  layer  to  the  cylinder  wall  is  equal  to  the  sensible 
heat  extracted  from  the  freezing  fluid,  together  with  the  latent  heat  of 


fusion  released  by  the  solidification  process  as  the  ice  layer  increases 
in  thickness.  The  cylinder  wall  is  assumed  to  have  negligible  thickness 
and  therefore  negligible  thermal  capacity  and  radial  resistance. 
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The  heat  balance  at  the  interface  may  thus  be  written 


k 


3.1-13 


l  3r 


6 
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3.2  Solution  for  Steady  State  Conditions 

When  a  superheated  liquid  is  flowing  through  the  freezing  zone, 
the  frozen  shell  of  ice  will  eventually  reach  a  steady  state  thickness 
for  every  axial  position  when,  at  the  interface,  the  heat  flux  in  the 
liquid  is  equal  to  the  heat  flux  in  the  ice.  The  resulting  steady  state 
interface  radius  will  be  used  as  a  reference  in  the  later  analysis  and 
hence  will  be  derived  separately.  It  will  be  assumed  that  the  temperature 
at  the  ice-water  interface  is  constant  and  equal  to  the  freezing  temperature. 

3.2-1  Temperature  Distribution 

The  approximate  form  of  the  energy  equation  in  the  ice  shell  is 


3. 2-1.1 


as  revealed  by  the  order  of  magnitude  analysis  in  Appendix  B.  The  appropri¬ 
ate  boundary  conditions  are 
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9T . 

Mtt-1)  =  hn(T  -T  ) 
i  7  O'  w  c' 

r  rQ 


W  =  Tf  • 


Normalizing  equation  3. 2-1.1  we  obtain 


9  <j>.  9(j) . 

- -  +  _! - L  =  n 

2  r*  9r*  u 
9r* 


where 


<f>-  = 


0, 

0 


)ci  ’  6i  Ti  "  Tf  ’  eic  Tf  " 


T. 


The  solution  of  the  reduced  equation  is  well  known  and  given  by 


d).j(r*)  =  c(z)  In  r*  +  D(z) 


3.2-1 .2 


The  normalized  boundary  conditions  are 


(^r*=l  =  •  Bic  [*i(1)  +  1] 


4>i  (5S*)  =  0 


where 


,r*  =  f  .6?  =  ^ 


c  k 


i 


0 


_s 

’s  r0 


If  4>.(1)  is  the  normalized  cylinder  wall  temperature,  then  the  solution 
for  the  temperature  in  the  ice  shell  is 
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^•(r*)  =  <(,.(1)  [1  -  . 


3.2-1  .3 


Using  equation  3. 2-1. 3  and  the  convective  boundary  condition  at  the  wall 
we  find  that 


Bi  In  6* 

V1'  =  1  -  Bi c  1  n6* 


3. 2-1. 4 


so  that  equation  3. 2-1. 3  becomes 


4>i  (r*) 


In  6*  -  In  r* 
s 

- In  6* 

Bi  s 


3. 2-1. 5 


3.2-2  Interface  Location 

The  energy  balance  at  the  ice-water  interface  under  steady  state 

d  6 

conditions  is  the  same  as  equation  3.1-13  except  that  ^  =  0  since  the 
interface  has  ceased  to  grow.  Thus 


9T 
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i  9r  <6. 


3.2-2. 1 


As  seen  from  equation  3.1-11, 
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l  3r  1 6 
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°°  c .  2 

Denoting  x  =  I  -J-R.'(l)  exp  (-A,  z*)  equation  3. 2-2.1  becomes 
j=0  L  J  J 
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-  k. 


9Ti.  2MVMx 


r'o  'f- 


i  9r  '6 


3. 2-2. 2 


Substituting  into  equation  3.2-2. 1  and  normalizing  we  obtain 
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2k0  Tn-Tf 

JL  .f.)  JL 
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3. 2-2. 3 


where  a  is  the  superheat  ratio. 

Substituting  equation  3. 2-1. 5  yields 


1n6s  =  BT7  'I 


3. 2-2-4 


2k 


where 


3  = 


£ 


aX 


Therefore 


5S  =  exp  (bT  -  • 

c 

It  is  worth  mentioning  that  for  the  Dirichlet  problem,  i.e. 

Bi  -►  °°  and  the  wall  temperature  constant,  this  result  reduces  to  the 
c 

result  given  by  Zerkle  and  Sunderland  [5],  as  expected. 
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3.3  Solution  for  Transient  Conditions 

3.3-1  Temperature  Distribution 

Consider  the  equation  for  the  heat  balance  at  the  ice-water 
interface  in  the  freezing  zone  as  given  by  equation  3.1-13  but  rewritten 
in  the  form 


2MW  x 


rU  9T . 

i  d  o  _  i  l 

6  ‘  piLi  dt  "  ‘  l  9r 


6* 


3. 3-1.1 


Considering  radial  displacement  normally  inward  from  the  wall,  i.e. 
r1  =  rg  -  r  and  6'  =  rQ  -  6 ,  then  the  above  equation  becomes 
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3. 3-1. 2 


Normalizing  equation  3. 3-1. 2  will  give 
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3. 3-1.3 


Further  rearrangement  yields 
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where 


Then 


2  2 
PiL.r  rn 
=  pip  0  1 

LkJe_Jt_  J 


" I  Cl  c 


de 

di 


3.3-1 .5 


For  a  growing  layer  of  ice  the  coefficient  of  the  right  hand 
term  in  equation  3. 3-1. 6  can  not  be  very  small  without  suppressing  the 
nature  of  the  problem,  i.e.  ice  formation.  Taking  the  coefficient  identi¬ 
cally  equal  to  unity  we  define  t  as 

2  2 

i  p.r  rn  c  . 

t  =  _L_  r.J—J3 _ 5 _ Eli 

x  Ste  L  k.  J 

l 

c  .0  . 

where  Ste  =  P.1  c1-  is  the  Stefan  number.  Lock  [11].  This  is  the  ratio 

Li 

of  the  sensible  heat  to  the  latent  heat.  The  normalized  interface 
equation  is  then 
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3. 3-1. 6 


Normalizing  the  conduction  equation  for  the  ice  shell  in  the  freezing 
zone,  as  given  by  equation  3.1-1,  in  a  similar  fashion  as  shown  previously 
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Substituting  for  tc,  the  normalized  conduction  equation  becomes 


82(b.  r  3cj),  8(f). 

— — L  _  - 8 - L  =  cfp  L 

8a2  1  -  V  9a  * 

The  normalized  initial  and  boundary  conditions  are  then 


3. 3-1. 7 


<f>.j(e.T)  =  0 
84, 

(ar)  rPBic  [*i(0>T)  +  1]  3-3-K8 

a — L) 

4>i  (1 »°°)  =  0 

where  <f>.(0,  t)  is  the  normalized  wall  temperature. 

For  the  problem  under  consideration  the  analysis  will  be  re¬ 
stricted  to  Ste  «  1,  i.e.  the  latent  heat  effects  dominates  over  the 
sensible  heat  effects.  For  an  ice-water  system  the  Stefan  number  is 
usually  less  than  unity  and  may  be  very  much  less.  For  this  reason,  the 
ice  interface  advances  very  slowly  inward  which  suggests  a  perturbation 
expansion  in  Ste  itself.  The  solution  of  equation  3. 3-1. 7  is  then  taken 
as 
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Substituting  into  equation  3. 3-1. 7  and  grouping  powers  of  Ste 
we  obtain  an  infinite  set  of  equations.  Here  we  will  only  consider  the 
asymptotic  solution  <f>.j0(a,T).  For  this  analysis,  we  will  call  the  temper¬ 
ature  a  "quasi-steady"  solution  since  the  steady  state  conduction  equation 
is  used.  The  temperature  distribution  is  thus  logarithmic  at  all  times. 

The  "quasi-steady"  temperature  distribution  is  given  by 


4>io(a,T)  =  -  In  (1  -  r  a)  +  B(t)  . 
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Using  the  boundary  conditions  (equations  3. 3-1. 8) 
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where  Q ( 0 ,x )  is  the  normalized  cylinder  wall  temperature  given  by 
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3.3-2  Interface  Location 

Using  the  zeroth  approximation  of  the  temperature  solution  in 
the  interface  equation  3. 3-1. 6  we  obtain 


de  _  9(N‘  o 
dx  3r 


1 


or  substituting  equation  3.3-1.12  gives 


3  .  3-2 . 1 
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p  c  p  p 

Further  rearrangement  yields 
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-  r  { - Si - }  .  3. 3-2. 3 
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P  Bic  p 

It  is  evident  that  the  solution  of  equation  3. 3-2. 3  will  predict  the 
interface  history  at  any  axial  position  for  a  given  flow  rate,  superheat 


ratio  and  Biot  number. 
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CHAPTER  IV 


SOLUTION  OF  THE  MERGED  PROBLEM 

4.1  Interface  Growth  and  Temperature  Distribution 

Numerical  results  were  obtained  for  the  theoretical  analysis 
using  an  IBM  360/67  computer,  in  the  Department  of  Computing  Science, 
University  of  Alberta.  Appendix  A  provides  the  program  which  determines 
the  eigenvalues,  eigenfunctions,  the  length  of  the  ice-free  zone,  the 
temperature  profile  at  the  entrance  to  the  freezing  zone  and  the  history 
of  the  interface  at  any  point  in  the  freezing  zone.  Axial  profiles  of 
ice  thickness  and  pressure  drop  were  also  calculated.  Representative 
results  are  shown  graphically  by  Figures  4.1  and  4.2. 

The  eigenvalues  of  equation  2.2-2  were  calculated  to  the  desired 
accuracy  by  substituting  an  evaluated  initial  guess  into  a  Newton-Raphson 
iterative  technique.  Equation  2.1-9  and  the  variable  x>  which  involve 
infinite  series,  were  evaluated  by  summing  the  first  forty  terms.  The 
integrals  in  equations  2.1-11  and  3.1-12  were  evaluated  by  means  of  the 
Simpson  rule,  while  the  differential  equations  for  determining  the  eigen¬ 
functions  and  interface  history  (equations  2.1-8  and  3. 3-2.3)  were  solved 
by  a  fourth-order  Runge-Kutta  method  of  forward  integration. 

In  the  analysis  it  was  assumed  that  the  properties  of  the  ice 
and  water  were  constant. 
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FIGURE  4.1  INTERFACE  POSITION  VS.  TIME 


=  5.15,  r0  =  0.125' 
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FIGURE  4.2  INTERFACE  POSITION  VS„  LENGTH  (FIXED  FLOW  RATE) 
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4.2  Pressure  Distribution 

The  axial  momentum  equation  within  the  liquid  in  the  freezing 
zone,  assuming  axi -symmetry  and  negligible  natural  convection  is  given 
by 
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•] 
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4.2-1 


where  p'  is  the  pressure  departure  from  the  hydrostatic  pressure  at  the 
wall,  i.e.  ||—  =  |^-+  pwg  .  The  correct  boundary  conditions  may  be 
stated 


vr(6,z)  =  vz(6,z)  =  0 

r  2 

v  (r,0)  =  2V[1  -  (f-)  ] 

0 

p'(r,0)  =  p0‘ 


From  an  order  of  magnitude  analysis  (Appendix  B)  of  the  radial  momentum 
equation  it  was  found  that  the  pressure  is  predominantly  a  function  of 
the  axial  position. 

Following  Zerkle  and  Sunderland  [5],  equation  4,2-1  may  be 
multiplied  by  r  and  integrated  from  r  =  0  to  r  =  6,  resulting  in 
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using  the  continuity  equation  and  integration  by  parts.  Substituting 
the  parabolic  velocity,  vz,  and  non-dimensional izing ,  the  above  equation 
becomes 


where 


dp*  _  16  d6*  .  32Pr  ri  ,  ,d6,2n 

dz*  -  -  36*5  dz*  6*4  [1  W  ] 
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then  upon  integration  this  becomes 


p*(z*) 


— —r  (1  -  6*4)  +  32Pr  \  ^ 

36*  0  6* 


4.2-3 


The  quantity  6*  in  this  expression  is  obtained  from  previous  results,  i.e. 

<5*  =  1  -  r  e,  where  e  is  calculated  for  a  specific  time  and  tube  radius. 

P 

A  representative  result  is  shown  graphically  in  Figure  4.3. 
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FIGURE  4.3  PRESSURE  DROP  VS  LENGTH  (FIXED  FLOW  RATE) 


PART  II 

EXPERIMENTAL  INVESTIGATION 


CHAPTER  V 


DESIGN  OF  APPARATUS 

The  experimental  apparatus  was  considered  as  two  systems;  one 
for  controlling  the  external  boundary  condition  and  the  other  for  the 
internal  fluid  flow  (see  Figure  5.1).  The  system  controlling  the  external 
boundary  condition  consisted  primarily  of  the  reservoir  and  piping  for 
the  transport  of  the  cooling  agent  (a  mixture  of  50%  water  and  50% 
ethylene  glycol)  to  the  test  section  cooling  jacket.  The  coolant  was 
stored  in  a  closed  reservoir,  1  -  1/2  feet  in  diameter  and  2  feet  high, 
which  contained  a  cooling  coil  connected  to  one  of  two  vapor-compression 
refrigeration  units.  This  particular  unit  had  a  capacity  of  one  ton  per 
hour  enabling  the  coolant  to  be  maintained  at  a  constant  temperature  of 
-15°F.  The  reservoir  was  connected  to  the  cooling  jacket  through  a  one 
inch  diameter  copper  tubing  recircul ati ng  system,  containing  a  motorized 
three-way  valve  and  pump. 

A  resistance  bulb  thermometer  (L  7033A  Balco),  which  was  located 
about  one  foot  downstream  of  the  three-way  valve,  was  used  in  conjunction 
with  a  resistance  thermometer  controller  (Honeywell  R7087D)  which  con¬ 
trolled  the  motor  (Honeywell  M904  Modutrol  Motor)  linked  to  the  three-way 
valve  (Honeywell  Series  1616  three-way  valve).  The  system  pump  recircu¬ 
lated  the  coolant  from  the  jacket  either  into  the  main  reservoir  or  to 
the  three-way  valve  which  regulated  the  mixing  of  the  -  15°F  coolant 
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1,2  THREE-WAY  VALVES 


3,4,5, 6  GATE  VALVES 
7  FLOW  METER 


FIGURE  5.1  SCHEMATIC  OF  EXPERIMENTAL  APPARATUS 
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from  the  reservoir  in  order  to  maintain  a  constant  temperature  set  on 
the  resistance  thermometer  controller.  The  reservoir  and  piping  were 
insulated  to  prevent  heat  leakage  into  the  system. 

The  second  system,  used  for  the  supply  and  temperature  control 
of  the  water,  consisted  of  two  constant  head  tanks,  supply  and  discharge 
piping,  motorized  three-way  valve,  flow  meter,  entry  length  and  test 
section . 

The  water  in  the  larger  header  tank  is  cooled  to  the  desired 
temperature  by  a  coil  around  which  an  ice  bank  of  fixed  thickness  is 
formed.  The  ice  bank  was  used  to  maintain  the  water  at  a  given  temperature 
below  the  water  main  temperature.  The  coil  was  connected  to  the  second 
refrigeration  unit  having  a  capacity  of  one  half  ton  per  hour.  The  water 
in  the  second  and  smaller  header  tank  was  maintained  at  the  water  main 
temperature.  Water  from  both  header  tanks  is  then  piped  to  the  motorized 
three-way  valve  where  the  "cold"  water  is  trimmed  to  the  appropriate 
temperature  by  mixing  with  the  "warm"  water  from  the  small  tank.  From 
the  mixing  valve  the  water  is  fed  through  1  inch  diameter  copper  tubing 
to  a  level  of  eleven  feet  below  the  water  level  in  the  tanks,  passing 
through  a  flow  meter  (Fisher  and  Porter  Ratosight)  and  a  flow  controlling 
gate  valve.  The  water  was  then  channeled  into  a  six  foot  section  of 
three  inch  diameter  copper  tubing,  the  last  three  feet  of  which  served 
as  the  test  section:  the  latter  was  surrounded  by  the  previously  mentioned 
cooling  jacket.  The  water  could  be  discharged  from  either  the  test  section 
or  the  flowmeter.  The  large  header  tank  and  all  the  piping  were  insulated. 


J  i<  i  '  -'ns  t*i  i'Pv;  »"  ,  'i'  n  (' "  ba  .  •  •  'Cm  .  •  (i  qrq 


1  .•  '2 »r  i  beioo.  >  i  i  '  b.n  .-»j  :  f  d  »i  v  ortl 


ar  a<>9n>tofdd  b9X  rd  do  jJnsd  sor  ns  rtorriw  onuons  ffoo  s  9*iudfi'T9qfn9d 

9nudf>'  >  d  r.9vle  6  is  isi&vi  9dd  nrsdo  sr.i  oi  b9au  t  i  ined  9or  9rtT  .bgrrnod 
nrid  oj  b9+D9nno3  asw  C roo  <  .enudt<*tsq  ud  r -gr>  ngdfiw  9dd  wof9d 


j  j  nr  6n«  >vfNr  bf  i  n  :  nw  vx  ->  u  1  't 


o*  b?  ; ■ q  narli  :r  a^lnsj  -\sissfi  ridod  movr  nsd&W 
9d6 rnqonqqfi  arid  od  bominr'i:  ar  n9d6w  “broa"  9dd  on9dw  9vfgv  \qfiw-99ndd 


gnraasq  ,ann,d  add  nr  fnvof  n9dsw  9dd  woFsd  d99d  ns/9f9  do  f9V9f  6  od 


o  991  99nnd  d2£  and  ,q  rdji  i9qooo  n9d9mfirb  donf  99idd 


. do^iojs t  enf-fooo 


CHAPTER  VI 


INSTRUMENTATION 

The  prime  instrumentation  problem  involved  the  measuring  of 
the  ice  thickness  inside  a  long  vertical  pipe  surrounded  by  a  completely 
insulated  cooling  jacket.  Gort  [8],  met  this  difficulty  with  a  stain¬ 
less  steel  probe  with  mechanical  arms  to  detect  the  ice-water  interface. 

A  potentially  superior  method  demonstrated  by  Bailey  and  Dula  [12]  employs 
ultrasonics  and  it  was  decided  to  adopt  this  approach.  The  ultrasonic 
equipment  (Branson  Instruments  Sonoray  600)  was  operated  as  a  pulse-echo 
device.  A  single  1/2  inch  cubical  immersion  transducer  facing  in  a 
radial  direction,  attached  to  a  1/2  inch  O.D.  stainless  steel  probe, 
acted  as  both  transmitter  and  receiver.  The  probe  was  located  along  the 
axis  of  the  test  section,  with  the  time  lag  between  the  echo  from  the 
tube  wall  and  the  echo  from  the  ice-water  interface  determining  the  ice 
thickness.  An  analog  signal  of  the  ice  thickness  was  obtained  through 
the  use  of  an  electronic  timing  gate  (Branson  Time  Analog  Gate):  the  time 
lag  represented  precal ibrated  ice  thickness.  The  calibration  was  ac¬ 
complished  using  equivalent  ice  thicknesses  in  aluminium  and  then  checking 
and  recalibrating  the  gate  using  a  previously-used  apparatus  for  studying 
planar  ice  growth  (see  Gunderson  [13]).  An  illustration  of  the  pulse- 
echo  signal  obtained  is  shown  in  Figure  6.1.  The  accuracy  of  the  instrument 
is  0.005  inches  or  better  depending  on  the  position  of  the  interface.  Since 
the  ice  was  not  isothermal,  thus  modifying  the  velocity  of  sound,  there  is  an 
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FIGURE  6.1 


PULSE  -  ECHO  OPERATION 
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error  inherent  in  the  readings.  This  error  was  found  to  be  very  small 
since  the  temperature  gradients  were  approximately  the  same  during  the 
experimental  runs  as  those  during  the  calibration  procedure. 

The  probe  was  able  to  move  up  and  down  in  the  test  section  by 
means  of  a  reversing  D.C.  motor  driving  a  lead  screw,  to  which  the  one- 
half  inch  diameter  probe  was  attached.  This  position  of  the  probe  was 
then  monitored  by  a  forty  turn  potentiometer:  a  change  in  potential  from 
a  prescribed  datum  indicated  the  distance  of  the  probe  from  the  test 
section  entry.  A  signal  from  the  potentiometer  along  with  the  analog 
signal  from  the  ultrasonic  gate  was  then  used  to  drive  an  X-Y  recorder, 
with  the  ice  thickness  plotted  on  the  vertical  axis  and  the  axial  position 
on  the  horizontal  axis.  By  this  means  it  was  hoped  that  an  axial  profile 
of  the  ice  thickness  at  any  time  could  be  generated. 

In  order  to  maintain  the  desired  coolant  temperature,  a  re¬ 
sistance  bulb  thermometer  was  located  in  the  coolant  recirculating  line, 
as  outlined  previously.  This  resistance  bulb  thermometer  was  connected 
to  the  temperature  control  mechanism  which  compared  the  thermometer's 
signal  to  a  preset  control  valve.  The  error  signal,  if  any,  energized 
the  motor  controlling  the  valve  movement  thus  opening  the  valve  for 
either  the  -  15°F  coolant  or  the  warmer  coolant  recirculated  from  the 
jacket,  depending  on  the  polarity  of  the  signal.  The  three-way  valve 
operation  for  control  of  the  water  temperature  was  similar,  except  that 
an  iron-constantine  thermocouple  sensor  was  used  and  the  valve  movement 
was  controlled  manually  by  the  use  of  a  three-way  toggle  switch. 
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Further  instrumentation  consisted  of  a  series  of  iron-con- 
stantine  thermocouples,  with  a  resolution  of  ±  0.25°F,  located  on  the 
outside  wall  of  the  test  section  at  discrete  positions  and  two  thermo¬ 
couples  located  in  the  center  of  the  tube  six  and  nine  inches  before  the 
test  section.  Readings  were  obtained  during  the  test  period  using  a 
low  sensitivity  sixteen  channel  recorder  (Brown-Honeywel 1 ) ,  with  a  range 
of  -  50°F  to  +  100°F. 
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CHAPTER  VII 


EXPERIMENTS 


7.1  No-Flow  Tests 

A  series  of  tests  were  conducted  without  water  flowing  through 
the  test  section  in  order  to  determine  the  reproducibility  of  the  results 
and  make  a  comparison  with  theory  and  other  experiments  (see  Gort  [8]). 

Before  each  series  of  tests  the  water  in  the  large  header 
tank  was  cooled  over  a  prolonged  period.  The  test  section  was  then 
completely  drained  and  the  cooling  jacket  recirculating  pump  started 
with  the  coolant  temperature  control  set  at  the  desired  temperature. 

After  waiting  to  ensure  the  stabilization  of  the  temperature  of  the 
coolant  in  the  cooling  jacket, the  valve  to  the  test  section  was  held 
open  until  water  was  seen  exiting  from  the  discharge  line  of  the  test 
section.  The  valve  was  then  closed  and  the  time  noted.  This  time  was 
now  considered  to  be  the  datum,  i.e.  t  =  t  =  0.  The  ultrasonic  timing 
gate  was  then  adjusted  carefully  along  with  the  gain  on  the  ultrasonic 
pulse  generator  so  as  to  trigger  on  the  ice-water  interface.  Since 
the  analog  gate  output  for  one  inch  of  ice  was  1 0V ,  readings  of  ice 
thickness  were  obtained  either  from  the  X-Y  recorder  or  a  digital  volt¬ 
meter.  The  wall  temperatures  were  then  taken  at  approximately  the 
same  time.  This  procedure  was  continued  until  the  ice-water  interface 
was  within  a  fifth  of  an  inch  of  the  transducer  face.  The  coolant  temperature 
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control  was  then  set  past  the  freezing  point  and  warm  water  passed  through 
the  test  section  thus  melting  the  ice.  After  the  ice  was  completely 
melted  the  warm  water  was  shut  off  and  the  test  section  drained.  The 
coolant  temperature  control  was  reset  to  the  desired  temperature  in 
preparation  for  the  next  experiment. 

7.2  Laminar  Flow  Tests 

The  next  series  of  tests  was  conducted  to  determine  the  in¬ 
fluence  of  flow  and  superheat  on  the  rate  of  ice  formation  and  make  a 
comparison  with  theoretical  predictions. 

With  the  coolant  temperature  control  set  past  the  freezing 
point,  the  flow  was  then  directed  into  the  test  section  where  its  flow 
was  measured  and  the  amount  of  superheat  set  using  the  entry  region 
thermocouples.  Once  the  flow  and  superheat  was  established  and  fixed, 
the  water  was  then  discharged  from  the  flowmeter.  Then  the  procedure 
for  preparing  the  test  section  was  executed,  as  in  the  no-flow  tests. 

The  desired  temperature  of  the  coolant  was  set  on  the  controller,  the 
pump  started  and  allowance  made  for  the  coolant  temperature  into  the 
jacket  to  reach  a  steady  value.  With  the  system  in  a  steady  condition, 
the  flow  was  then  switched  back  to  the  test  section  with  the  datum  time 
noted.  Thickness  and  temperature  readings  were  then  taken  with  the  flow 
being  checked  periodically  by  weighing  the  discharged  water  over  a  definite 
period  of  time.  The  end  of  each  test  was  determined  by  the  prevailing 
flow  and  superheat  ratios.  Each  test  was  terminated  when  there  was  no 
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significant  change  in  ice  thickness  after  fifteen  minutes  or  the  ice 
grew  within  a  fifth  of  an  inch  of  the  transducer  face:  a  test  was  termi¬ 
nated  if  the  superheat  of  the  water  rose  substantially  making  the  sub¬ 
sequent  results  useless  for  comparison.  The  test  section  was  then 
prepared  for  the  next  experiment  in  the  same  manner  as  the  no-flow 
tests:  coolant  temperature  control  set  past  the  freezing  point,  warm 
water  used  to  completely  melt  the  ice,  and  the  test  section  drained. 
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DISCUSSION  AND  CONCLUDING  REMARKS 
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CHAPTER  VIII 


DISCUSSION 


8.1  No-Flow  Problem 

This  series  of  experiments  was  conducted  without  flow,  as  dis¬ 
cussed  previously,  to  determine  the  reproducibility  of  the  results  and 
to  make  a  comparison  with  other  experiments  and  theory.  The  theory  for 
the  no-flow  problem  in  the  absence  of  superheat  is  essentially  the  same 
as  for  the  problem  where  a  liquid  is  flowing  with  its  bulk  mean  temperature 
at  the  freezing  point. 

A  series  of  four  tests  were  run,  with  the  water  temperature 
maintained  as  close  to  the  freezing  temperature  as  possible.  The  results 
were  compared  with  those  obtained  by  Gort  [8],  who  assumed  the  wall 
temperature  difference  decreased  according  to  a  power  law.  Using  Gort's 
approach,  i.e.  a  Dirichlet  boundary  condition,  the  results  obtained  for 
the  no-flow  tests  were  in  close  agreement  (see  Figure  8.1)  and  reproducible. 

All  readings  were  taken  at  twenty  seven  inches  from  the  leading 
edge  of  the  test  section.  This  was  done  since  the  ultrasonic  transducer 
was  then  opposite  two  thermocouples  which  were  used  to  determine  a  mean 
heat  transfer  coefficient  at  this  point.  During  one  of  the  tests  a 
limited  axial  traverse  of  the  interface  was  undertaken  using  the  digital 
voltmeter:  that  is,  four  inches  either  side  of  this  position.  It  was 
found  that  the  profile  did  not  change  appreciably  (less  than  0.008  inches). 


49 


ni'  JUraHCi 


H0I22l'32ia 


N 


229  )  \  d6b9'tqq6  9gn5fl3  don  brb  sfr^oiq  9dd  daqd  bnuo^ 

- 


50 


oo\o 


FIGURE  8.1  NO  -  FLOW  RESULTS 
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This  position  was  then  considered  to  be  a  representative  location  for  the 
theory  since  the  rate  of  growth  is  independant  of  axial  position  for 
this  particular  case.  A  mean  heat  transfer  coefficient  was  then  calcu¬ 
lated,  using  the  wall  thermocouple  readings  and  the  corresponding  coolant 

temperature,  from  equation  3.3-1.11.  This  revealed  that  Bi  =  5.15,  27 

c 

inches  from  the  leading  edge. 

The  theoretical  curve  of  6/r^  versus  t  is  plotted,  along  with 
the  experimental  results,  for  the  evaluated  convective  boundary  condition. 
The  comparison  given  in  Figure  8.2  shows  that  there  are  three  possible 
factors  which  would  be  responsible  for  the  discrepancy  between  the  experi¬ 
mental  data  and  the  corresponding  theory.  The  primary  factor  would  be 
the  superheat  present  in  the  water.  The  greater  the  superheat,  the  greater 
the  deviation  to  be  expected  as  the  superheat  must  first  be  removed  from 
the  water  thus  causing  the  experimental  curve  to  be  beneath  the  theoretical 
curve.  The  second  factor  would  be  the  non-zero  Stefan  number.  Since 
the  Stefan  number  is  non-zero,  the  rate  of  ice  growth  would  be  slower 
since  it  is  this  non-dimensional  group  which  reflects  the  significance 
of  the  sensible  heat  of  the  ice.  This  effect  is  small  as  noted  in 
Chapter  9.  Finally  there  is  the  consideration  of  natural  convection 
inherent  in  the  test  section  due  to  the  initial  superheat  present  in  the 
water.  At  the  ice-water  interface  cooling  is  continually  taking  place 
thereby  lowering  the  temperature  of  the  water  toward  the  freezing 
temperature,  and  thus  decreasing  the  value  of  its  density.  Since  the 
water  was  always  below  its  inversion  point  (approximately  40°F),  a  circu- 
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lation  was  then  initiated  with  the  colder  water  ascending  along  the 
interface,  and  the  warmer  water  descending  at  the  center  of  the  pipe. 

The  inclusion  of  natural  convection  would  tend  to  decrease  the  rate  of 
growth  of  the  ice,  thus  again  causing  the  experimental  points  to  locate 
below  the  theoretical  curve.  Since  no  experiments  were  carried  out  to 
determine  the  magnitude  of  this  effect,  an  estimate  of  this  discrepancy 
compared  to  the  actual  observations  was  not  possible.  The  water  temperature 
was  always  brought  as  close  to  32°F  as  possible  (e.g.  34°  -  36°F)  thus 
making  the  effect  of  natural  convection  as  small  as  the  system  would 
allow. 

From  superimposing  Figures  8.1  and  8.2  it  is  evident  that  the 
theoretical  predictions  from  different  viewpoints  (i.e.  Dirichlet  and 
convective  boundary  conditions)  of  the  same  problem  are  in  close  agreement. 
This  reveals  that  the  interface  history  can  be  predicted  accurately  using 
either  theory  for  the  no-flow  tests. 

8.2  Laminar  Flow  Problem 

A  series  of  tests  were  run  for  different  superheat  ratios  and 
Reynolds  numbers  with  the  results  being  compared  with  those  predicted  by 
the  theory  developed  earlier. 

The  results  are  presented  in  Figures  8.3  -  8.7:  the  Reynolds 
number  in  each  test  was  based  on  the  pipe  diameter  at  the  entrance  to  the 
test  section.  Originally,  axial  profiles  were  to  be  taken  at  designated 
time  intervals  during  a  test.  After  considerable  testing,  it  was  found 
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FIGURE  8.7  EXPERIMENT  WITH  FLOW  Re  =  345 
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that  the  ultrasonic  equipment  was  unable  to  pick  up  a  continuous  return 
of  echoes  necessary  to  operate  the  electronic  gate  required  to  obtain  a 
readable  plot.  The  digital  voltmeter  was  then  used  to  obtain  the  necessary 
readings  at  the  selected  position.  As  with  the  no-flow  tests,  the  trans¬ 
ducer  was  placed  twenty-seven  inches  from  the  leading  edge  for  two 
reasons.  The  first  one  is  the  desirability  of  being  opposite  the  thermo¬ 
couples  to  facilitate  the  checking  of  the  mean  heat  transfer  coefficient. 

The  second  reason  is  that  the  method  of  discharge  of  the  water  from  the 
test  section  to  the  drain  caused  the  least  disturbance  at  this  location. 

This  was  demonstrated  by  a  limited  number  of  axial  and  azimuthal  traverses 
of  the  probe  done  in  preliminary  testing.  The  azimuthal  traverses  also 
demonstrated  that  the  assumption  of  axi-symmetry  was  valid  since  the 
variations  in  ice  thickness  were  only  ±  0.008  inches.  It  should  be  noted 
that  the  ice  was  still  growing  when  the  traverses  were  made  so  the  vari¬ 
ations  were  probably  smaller.  During  one  preliminary  test  it  was  found 
that  an  ice-free  zone  did  not  exist  in  the  test  section  as  expected.  Since 
axial  profiles  were  not  obtainable,  a  check  to  determine  if  an  ice-free 
zone  existed  was  not  possible.  Since  the  system  Nusselt  number  was  large 
(20.79)  it  was  assumed  that  an  ice-free  zone  did  not  exist  over  the 
temperature  and  flow  ranges  tested  (see  Figure  2.2). 

There  are  several  factors  which  would  explain  why  the  experi¬ 
mental  points  deviated  from  the  theory.  First  of  all,  in  Figures  8.3 
and  8.4,  the  flow  was  high  enough  that  the  local  Reynolds  number  was 
either  initially  above  the  critical  Reynolds  number,  taken  as  Recrit  ~  2000, 
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or  increased,  due  to  increasing  ice  thickness,  to  a  point  where  it  passed 
the  critical  Reynolds  number.  Since  the  flow  is  then  no  longer  laminar, 
the  theory  would  no  longer  be  applicable  beyond  this  ice  thickness.  Since 
turbulent  flow  results  in  an  increased  heat  transfer  rate,  it  is  reason¬ 
able  to  assume  that  the  steady-state  ice  thickness  would  be  much  lower 
than  that  predicted  by  the  laminar  flow  theory.  For  these  higher  flow 
rates,  this  increase  in  local  Reynolds  number  to  a  point  exceeding  the 
critical  value  would  be  the  dominant  factor  in  the  discrepancy  between 
the  theory  and  experiment.  The  experimental  curves  of  Figure  8,4  are 
in  close  agreement  with  those  predicted  by  theory  in  that  region  where 
the  theory  is  applicable. 

Next  consider  the  situation  where  the  initial  Reynolds  number 
was  low  enough  to  preclude  turbulent  flow.  Here  the  discrepancies  between 
theory  and  experiment,  shown  in  Figures  8.5  to  8.7,  appears  to  be  a  result 
of  the  neglect  of  the  effect  of  natural  convection  at  these  lower  flow  rates 
This  possibility  is  supported  by  an  order  of  magnitude  analysis  of  the 
momentum  equation  shown  in  the  Appendix.  The  circulation  now  initiated 
is  quite  different  from  that  in  the  no-flow  case,  as  the  water  temperature 
is  now  above  the  inversion  point.  The  natural  convection  current  alone, 
without  a  superimposed  forced  flow,  would  then  move  upward  at  the  inter¬ 
face  and  at  the  center  of  the  pipe,  though  downward  in  a  region  between. 

This  effect  is  quite  complicated  and  it  is  difficult  to  say  precisely 
what  effect  it  did  have  on  the  flow  tests.  We  know  the  effect  of  natural 
convection  will  diminish  as  the  water  flows  through  the  test  section 
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since  its  mean  temperature  approaches  the  freezing  temperature.  There¬ 
fore  the  effect  of  natural  convection  would  be  the  greatest  in  the  thermal 
entrance  region.  This  appearance  of  natural  convection  causes  heat 
transfer  to  be  greater  than  theoretical  predictions  with  the  result  that 
the  radius  of  the  ice-water  interface  is  also  greater  than  the  theoretical 
predictions.  Also,  as  the  ice  thickness  increases  along  the  tube,  due 
to  the  mean  temperature  of  the  water  approaching  the  freezing  temperature, 
the  flow  area  is  constricted  resulting  in  the  water  accelerating  and 
thereby  reducing  the  effect  of  natural  convection  in  comparison  with  forced 
convection.  Since  no  tests  were  conducted  to  determine  the  effect  of 
natural  convection  it  is  extremely  difficult  to  say  what  quantitative 
discrepancies  are  involved,  except  through  the  order  of  magnitude  analysis. 

A  comparison  was  made  for  one  of  the  flow  tests  with  a  similar 
experimental  test  made  by  Gort  [8],  shown  in  Figure  8.5.  As  expected, 
the  data  obtained  by  Gort  revealed  that  the  ice  grew  faster  initially  and 
then  grew  slower  as  the  superheat  ratio  increased,  as  opposed  to  the 
experimental  test  run  at  a  constant  superheat  ratio. 

For  all  set  of  curves,  the  question  could  arise  as  to  whether 
the  experiments  were  outside  the  hydrodynamic  entry  length  where  the 
velocity  profile  is  fully  developed  and  the  theory  applicable.  As  shown 
by  the  theoretical  interface  profiles  at  prescribed  times,  the  hydro- 
dynamic  entry  length  increases  as  the  ice  thickness  increases,  i.e.  the 

hydrodynamic  entry  length  is  a  function  of  time.  This  can  be  demonstrated 

d  (5 

by  Figure  4.2,  the  position  of  a  particular  interface  slope,  ,  moves 
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downstream  as  the  time  progressed.  If  the  point  at  which  the  experi¬ 
mental  results  were  taken  was  to  be  in  the  entry  region,  or  the  entry 
region  increased  in  time  to  include  this  location,  then  an  appreciable 
error  in  the  theoretical  predictions  would  arise. 

For  all  experimental  cases,  the  effect  of  the  non-zero  Stefan 
number  would  also  have  to  be  considered.  This  would  decrease  the  rate 
of  growth  of  the  ice  due  to  the  inclusion  of  the  sensible  heat  of  the 
system.  The  neglect  of  the  unsteady  term  in  the  energy  equation  of  the 
liquid  has  a  very  small  effect  on  the  theoretical  predictions  for  these 
experiments,  as  noted  in  the  Appendix,  but  would  be  important  for  tests 
involving  much  larger  superheat  ratios. 

Noting  that  the  average  Biot  number  for  the  experimental 
apparatus,  at  the  point  considered,  was  established  using  results  from 
the  no-flow  tests.  Deviations  from  the  mean  value  could  be  expected  in 
some  experiments  due  to  the  coolant  temperature  not  remaining  constant. 

This  would  fluctuate  the  experimental  points  depending  on  the  temperature 
variation.  Not  only  will  the  coolant  temperature  fluctuate  entering  the 
test  section  but  it  could  vary  with  axial  position.  In  a  counter-flow 
heat  exchanger,  as  is  the  experimental  test  section,  the  effect  would 
be  to  increase  the  coolant  temperature  as  it  approaches  the  leading  edge 
of  the  test  section.  This  would  decrease  the  predicted  steady  state 
radius  as  well  as  the  complete  history  of  the  interface.  This  effect 
would  be  propagated  downstream  to  affect  results  everywhere.  The  deviations 
caused  were  difficult  to  assess  since  axial  profiles  were  unobtainable, 
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but  would  likely  be  noticeable  at  higher  superheat  ratios  and  flow  rates. 
This  is  possibly  substantiated  in  comparing  Figures  8.5  and  8.7. 

In  the  theoretical  solution  of  the  ice-free  zone,  whenever  z+ 

e 

is  small,  i.e.  z^  <  0.001,  a  large  number  of  terms  in  the  series, 
equation  2.1-10,  must  be  taken.  Since  the  computer  program  for  the  ice- 
free  zone  did  not  facilitate  a  great  many  terms  it  was  assumed  that, 
for  these  conditions,  z  +  =  0  and  a  block  temperature  profile  was  then 

c 

used  at  inlet  to  the  freezing  zone. 

For  the  plane  problem  Savino  and  Siegel  [7]  observed  that  the 
experimental  transient  ice  thicknesses  were  within  15  percent  of  the 
theoretically  predicted  values.  The  experimental  values  given  here 
(Figures  8.1  to  8.7)  are  generally  within  this  departure.  Considering 
the  experimental  curves  overall,  it  has  been  demonstrated  that  the  larger 
the  superheat  ratio  present  in  the  flowing  water  the  more  it  retards  the 
growth  of  the  ice  resulting  in  a  smaller  steady  state  radius.  Since  the 
Reynolds  number  is  an  indication  of  the  velocity  of  the  water  flowing 
through  the  test  section,  it  is  shown  that  the  increased  velocity  of  the 
water  results  in  a  slower  rate  of  growth  due  to  more  superheat  per  unit 
time  being  carried  over  into  the  test  section.  Both  these  results  confirm 
expectations . 
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CHAPTER  IX 


CONCLUDING  REMARKS 


9.1  Conclusion 

This  thesis  has  presented  a  theoretical  and  experimental 
analysis  of  liquid  solidification  in  a  vertical  tube,  with  prescribed 
convection  as  the  external  boundary  and  interface  condition,  in  the 
laminar  flow  regime. 

The  theoretical  approach  considered  only  the  asymptotic  solu¬ 
tion  for  the  temperature  profile  in  the  growing  ice  phase  since  for  ice- 
water  systems  the  Stefan  number  is  usually  much  less  than  unity.  The 
use  of  this  temperature  profile  in  the  interface  equation  yielded  an 
asymptotic  solution  for  the  history  of  the  interface.  The  effect  of 
neglecting  the  sensible  heat  effects  in  the  ice,  has  been  shown  by  Lock  [11] 
to  be  small  for  the  plane  problem.  Using  the  experimental  Stefan  numbers 
and  superheat  ratios  from  the  cylindrical  tests,  it  is  not  unreasonable 
to  expect  that  approximately  the  same  error  would  be  introduced  in  the 
cylindrical  system,  at  least  during  the  early  stages  of  ice  formation 
when  the  forming  phase  would  be  nearly  planar. 

The  statement  of  overall  heat  balance  at  the  interface  was 

found  to  be  dependent  upon  two  dimensionless  groups  contained  in  the 

ko 

variable  3:  the  superheat  ratio  a,  and  the  conductivity  ratio  .  For 
a  given  liquid  (i.e.  k£  and  k.  known),  the  greater  the  value  of  a  ,  the 
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greater  the  rate  of  freezing  and  thicker  is  the  ice  at  a  given  axial 
position  z. 

The  majority  of  the  theoretical  predictions  were  well  borne 
out  by  the  experimental  results  obtained  within  a  region  where  the 
theory  was  applicable.  A  comparison  of  theoretical  and  experimental 
results  of  —  versus  t,  for  all  experiments,  generally  revealed  agree¬ 
ment  to  within  20%  or  better.  Discrepancies,  although  never  large, 
were  greatest  as  the  ice-water  interface  approached  the  center  of  the 
tube.  On  occasions,  the  signals  received  by  the  ultrasonic  transducer 
became  random  to  a  degree,  e.g.  ±  0.020  inches,  probably  due  to  the 
increasing  concavity  of  the  interface  resulting  in  the  reception  of 
multiple  echoes.  As  mentioned  previously  the  effect  of  natural  convection, 
as  well  as  the  effect  of  being  in  the  hydrodynamic  entry  region  where 
the  inertia  terms  are  of  importance,  would  have  affected  the  results. 

A  considerable  amount  of  work  was  done  on  the  experimental 
apparatus  to  try  to  obtain  axial  profiles  which  are  quite  important  in 
making  a  complete  discussion  of  the  theory  and  experiment.  The  various 
modifications  to  the  apparatus  did  produce  (or  lead  to)  limited  axial 
profiles  using  the  X-Y  plotter.  A  curve  is  shown  in  Figure  9.1  with  the 
results  unplottable  at  a  distance  of  24  inches  from  the  leading  edge, 
as  the  illustration  shows. 

9.2  Suggestions  for  Further  Work 

The  analysis  is  limited  to  the  history  and  profile  of  the  ice 
growth  in  a  vertical  tube  neglecting  natural  convection,  inertia  terms 
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FIGURE  9.1  AXIAL  PROFILE  TRIAL 


67 


and  the  unsteady  energy  term  in  the  fluid.  A  logical  extension  would 
be  to  include  these  effects  in  the  analysis  as  well  as  the  effect  of 
turbulent  conditions.  Also,  other  positions  such  as  a  horizontal  or  an 
inclined  tube  could  be  studied  in  conjunction  with  natural  convection 
and  the  inertia  terms. 

It  was  also  difficult  to  make  any  comparison  in  the  experiments 

•  • 

reported  here  and  the  work  of  Ozisik  and  Mulligan  [14]  as  they  considered 
a  Dirichlet  boundary  condition,  i.e.  Bic  -*  ».  An  experimental  comparison 
could  be  made  if  the  apparatus  were  to  be  built  so  this  condition  could 
be  examined. 

The  ultrasonic  equipment  is  quite  suited  for  horizontal  and 
inclined  tube  investigations  as  the  transducer  does  not  have  to  be  located 
in  the  center  of  the  tube  to  obtain  accurate  readings.  The  experiments 
demonstrated  the  reproducibility  of  the  experimental  apparatus  and  the 
potential  of  the  ultrasonic  equipment. 
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FORTRAN  tv  PROGRAM  for  thf  MERGED  PROBLEM 


t  IN  Cnl*  4  TNOICATPS  THAT  THIS 

LINE  TS  A  PART  OF  THE  PREVIOUS  CARD 


C  CYLINDRICAL  ICE  FOR  MATT  ON 
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r  IFF  F-R.FF  7  r>N  f 

r 
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*  40) 

DIMENSION  X  ( 1  0?  )  ♦  Y  ( 1 0? )  *  V* (  1  OR  ) 

REAL  r'ATA(10R4) 
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F ( A) =  (  1 •O/BTPT )  +  ( ( CONST*SIN < ( A*PT/4#0 )-( ?.0**>I /3, 0) ) ) / 
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P  I  =  3 . 1  A 1  5P?6 

BinT=Bi*ricF/CH^n 

CON$T=  f  (  3.0**0«.  666667  )  ’fcGA'^MA  (  1. 337  333  ))/((  ?, 0**0. 37773 
*  73 ) *Q AMMA ( 0. 

7666^67’ ) ) 

T  =  0 

ncLw«o.oi 

WRTTr ( 106 ) 

VI P  T  T  F  f  6  t  5 1 1  ) 

VIP  T  TF  (  6, 307) 

WRTTF(6t50Q)  3TPT 
WP TTF (6,510)  TL 
VIP!  TF(  6, 500  ) 
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r  CALCULATE  I  MI  T I  At.  OOF  55  FOP  PHOT 
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A  =  0  F  L  W 
71  8*A+DFLW 

C !  J  M  7  —  A  B  5  (  F  {  A  )  +•  F  P  (  R  )  ) 
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2  0  A 1  =  A 
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CALI  PIN(RLAMO) 
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c  50LVF  T H F  DIFFERENTIAL  EQUATION  FOR 

r  j  mtpqp  at  F  history  bv  forward  integration 
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R  (  m,b  )=p  (  M,  X  )  +H/?.0 

o (Mf 4) =R ( M, i ) fH 
T(M)=0.0 

DR=  (  1 , 0-RP*R  (N,in 
WRTTR{6t10P)  R(N,1),T(N),°R 

X( N) -o ( Mf 1  ) 

Y  (  N  )  =  T  {  N  ) 
nn  i  w=i I  f 4 

XX  (  M)=1  .o*RP*R  (  NfNI) 

1  TP(N,M)=XX(M)*( GAM-ALOG(XX( M) H /f< 1 . O-BETA* (GAM- ALOG( X 

*  X(M)) ))*RP) 

OFLT  (N)s(H/6#0)*(  Tp(M,  n+(  ?.0*TP(N,  3)  )  M  ?.0*Tp(Nt3)  )+T 

*  P  f  N  ,  4  )  ) 

no  ?  f|-7  f  1  no 

R(N,  1  )=R(N-1  »n+H 
P(N,?)=R(N,l)+H/?.0 

Q  (NOI  =R  (N,  1  )  +H/7.0 
o  (  NJ,4)=R  (  isjt  i  )  *h 
T(  N) =T ( nj-1  ) ♦DELT ( M-1  ) 
nn  3  m~i  f  4 
X  X  (  M  )  -  1 . 0-RP*R (  N ,  M ) 
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7  Y(N!=T(N) 
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CALI  P*  OT  (0.0,0.0,73) 
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CALL  PLOT(6.0,l?.Ft?) 
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CALL  PLOT (i 7.0,0.0,33 ) 
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WRT  TC( 6, 1 05 ) 

WR I T  r ( 5 , 7  00) ( I  X  I  K  )  ,  Y  (  K  )  )  , K= 1  ,100) 

t  =o 


CALCULATE  AND  PLHT  AXIAL  PROFILES 
C  OF  R  ART  y  $  AND  PRESSURE 

r 

mop=mz*iqo 

RIGZ  =  FLOAT(Nf}P) *0.001 
CALL  PRFSS1 ( BTG7 , TIME ) 

WRTTc:(^t3fH  )  TTVF,Or| 

HO  A  M=MP,MTfy£ 

T  =  I  +  1 

77=eloaT( m) *0.00! 

X(T  )  -17 

Y( I )=DFLZ!ZZ,TIvf) 

WR!TF(6,300)  XU),Yf!) 

I Ff I .eq.LOO)  GO  to  s 

4  CONTINUE 

5  IF ( KONST. FO. 1 )  GO  to  oq? 

CALI  SCALF(XtS.O,)no,1 ,20.0) 

CALL  SCALE* Y ,6. 0 , 100 , 1 , 20. 0 ) 

CALL  AXIS (0.0,0. 0,70H0TMFNS I ONL FS^  LFNGTH,-?0, 8 . 0, 0.0, 

*  X( 101 ) ,X(102 
1 ),?0.Q) 

CALL  AXTS(0.0,0.0,20H  RADIUS  <R/po)  , ?o , 4. 0, q0. 0, 

*  Y ( 10) ) , Y( 10? 

1 1 , 70.0) 

CALL  LINE <X,Y,1 00 ,1,0,0) 

CALL  DR AW 7 (TIME, ZCR IT) 

CALL  PLOTfO. 0,0.0,?) 

WR  T  TF  (  6  ,  L  06  ) 

WRTTF(6,?00)((X(K),YtK))  ,K= 1,100) 

CALL  PLnT<0.0,14.0,?3 ) 

70?  WR  T  TF ( 6,302)  TImf,RO 

nq  6  1=1,100 

R  R =  v  (  T  ) 

W(T)=PRFSS(PR»RR*I) 

6  WRITC( 6,700)  X! t ) ,W( T ) 

I F ( KONST. FO. 1 )  GO  TO  PQO 
CALL  SCALF(W, 6.0, 100, ! ,70.0) 

C  A I  L  SCALE(X,8.0, 100, 1,20.0) 

CALI  AXIS  <0.0,0. 0,?0HD!NFNS IONL FSc  LENGTH, -70, 8 . 0, 0. 0, 

*  XI 101 ) ,X( 10° 

1 ) ,70.0) 

CALI  A X I S I  6. 0, 0 • 0 »  23HDIMFNS I PNL F SS  PRESSURE  ,73,6.0,90 

*  .0,WH01),V!( 

1 1 0  p )  , 70.0 ) 

TALL  I  INFIX, W, IOC*1 ,0,0) 

CALL  DRAW?! time, 7CRT T) 

CALL  PLnT  10. C,0.0,’ ) 

CALL  PLOT (l 2. 0,-1 4. 0,3) 

CALI  PI.  0T  (0.0,0.0,990) 
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CND  p j_ n T  ROUT INF 
WRITE!  6,106  ) 

WRTTp(6»000) ! ( X  !  K  )  ,  W  f  K  ) ) fK= 1,100) 

1.9  0  EHRM^t  (»  o  »  f  ]qx  ST^CE  ''ZCR!T««  IF  LpSS  than  0*001  THE 

*  0  Lnf  K  ppoft 

1 Lc  TL  =  *  *  F6. ? , *  ts  then  INPUT  INTO  TH^  TCr  PROBLFM* ) 
101  FORM  AT ( T  2  ) 

10?  FnRVAT(Rp^.3tpAs2, ?F6,o) 

104  FORMAT (♦ 0* , EX , 'CONSTANTS  KI  =  ' *  F4*  2  * ' *  KL  =  *,F6.3 

*  pp  =  »  f  f 

c6*of.f  mu  =  *  *F6,  3 ,  *  *  TF  *  *  *  F6 . ?  ,  *  ,  Tf_  =  'fF6.?,'f  tf 

*  =  '  *  F6*  ? ,  '  , 

f  RHH  =  •  t  Ffc.B/ ) 

105  E0RMAT(7F8*3,F1  2.6,c12.6) 

10^  FORM  AT ! ' 1 • ) 

107  FORM  AT { « 0» , ?X , »  ICE  INTERFACE  POSITION  AND  TI mf  ftp 

*  pc/nn  =  «  f 

?F9.6,!0Xt fPHHO/K  =  ' , F7. 0, lOXt *  Z*  =  *,F7.4) 

100  FORM  AT  f  '  ' , 1  OX , *  PFFTTION  RO-R/RO-RS  =  * * F6. 4 , 1  OX , * D T M 

*  ENS IONLESS 

7TTMC  =  *  ,fo.6t 1  OX, * R /RO  -  * , F6 • 4 ) 

700  FORM  AT ( *  * ,F6*4, Fl?*6  ) 

001  FORM  AT (»i  *  f  »  7  R /R0» ,1  OX,  * TIMp  -  «,FT*4,*  HRS 

*  ,  RO  =  * 

1 ,F6. 7, «  FT. » / 1 

700  FORMAT { *  1  * , •  7  PRESSURE* , Vox* 'TIME  =  '*F7.4*»  H 

*  R  S  ,  PO  = 

^  *  *  F  6*  3  *  •  FT../) 

407  FORMATf.  * ,10X, 12, 10X,F10*6,10X,F10*6) 

40  5  FORM  AT <  » i  » , i OX , *N» , 1  ox, *  FTGFNVALUF* , i 5X, »C ( M) * // ) 

40  A  format  ! *0' , 1  OX* 'N' *  1  OX, *-1 /?C(N) R  "  !1) *  //) 

50  0  FORMAT ( » 0* , 10X * ' T »  , 1  OX, *  F  TGENVALUF* // ) 

501  FORM  AT ( »  »  ,  10X* I?  ,10X ,F10. 6) 

F04  FORMAT  !»  O' *10X* • T ', 10X*  *CK ! t )» //) 

FO A  FORMAT! *  O' ) 

fo7  FORMAT (* 0* , 10X, » FRFF7 img  TEMPER ATURF  CHFC K  TF  =  *,F7. 

*  0 ) 

for  FORMAT! *  O' ,10X, • EIGFNVAUJFS  FOR  thf  TCF  frff  ZONF ' ) 

F0  °  FORMAT!*  *  * 10X* • FOR  »*RIOT»»  NO.  =  *,F8.3/) 

F 1  0  FORMat!  *0*  »10X»  *PLnFK  TMLFT  tFmpfraturf  PROFILE  =  *  ,  F6 

*  .V/) 

51.1  FORMAT!  «  0*  ,  BOX* 'CYLINDRICAL  T  F  c  FORMATION  WITH  CONVERT 

*  T vf  BOUNDARY 

1  CONDITIONS'//) 

*00  FORMAT !» O' , 10X* • ICE  WILL  NOT  FORM  UNTIL  * ' Z/RO*RF*PR *  * 

*  =  * , F8*  6 ) 

7Q 0  FORMAT !F7.4,?T4) 
ooo  FORM  AT ! 7 17* 7F7. 7 ) 

o 0 i  FORMAT!* O' ,?FX, *R/RD  vs.  (T-TF ) / ( TL-TF)  at  thf  I NL p 

*  T  TO  rMr  Tpc 
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FUNCTION  por^  i  PR  ,  op  f  T  ) 

COVMfN  MHTf  tf,  TCfTL,C!CF,CH20,8! 

C  H  M  M  H N  RHM(40)  ,RLAMW(40)  tCH<40)  tCC<40>  *CK(40)  ,RIR(  100 

*  i  ) 

COMMON  RTR  (  I  01  ,40)  ,  PRCS  (101  )fVM10’  ) 

CS-PPPS ( I +1 ) 

PSO=  po*pp 

riNM^o/^.oi^n  •n-R$q*p sen  /<r$q*rsq) 

PR  cS S -CON  +  "* ? 

RFTIJRN 

FNID 


M,"'-.  ''  *  4  ,T  "»r 

T  ,  .  ’  '  .  .  .  ’  .  T  ■ 

• ) '  ■  f  1  •  . '  ’  *1  '  .  i  .  . '  "  1 

< f 

{  t  \  ‘v  f  (  •  n  *>  '  «r  *  '  ^ 

t  '4  *  1  ‘  -“> 

»\(.  ?  '  "* 

rr-» 


FUNCTION  PEL  7 ( 7 , TQ ) 

HT  mcnjc  ton,  yx  (4)  ,TP(!  001,4)  ,  CELT  (1001  )  ,  R ( !  001 , 4 ) , T ( \ 00 1 

A  ) 

CnMMPNj  R  T  n  T  t  T  F ,  TC»TL  ,CTCE,CH?0,8T 

Common,  RL AM (40) ,RLAMW(40) ,CH{ 40 ) ,CC< 40)  ,C  K  (  40  )  ,  RL  R  (  1  00 
♦  1  ) 

COMMON  RTPnci,40)tPRBUOl),VM101) 
r  0 mm qm  CL  T ,RH0T ,  PP 
H=0 , 00 1 

CALI  PPF( Z, GAM, BETA, RS,RP) 

TE(R$.EQ.!.0)  on  TO  OR 
THETAC=TF-TC 

TCOMST-=  (  CL  T  /THET  AC  )  *R.  H0T*RP*RP*Rn*RO/r  t  CF 
T  z  =  T  q  /  t  C  0  N  S  T 
N  =  ! 

R  (  N ,  !  )  =  0 . 0 
P (N,?)=P (N,l ) +H/2.0 
R<NV?»=R(N,  1)  +H/2.0 
R  (  N  ,  4  )  =R  (  M ,  1  )  +H 
T{ N)=0.0 
TP^T=T(M) 

TF( TFET.GF.T7)  GO  TH  20 
nn  i  m=| f 4 

KKf M1=1.0-RP*R(NfMJ 

1  TP(m,M)=xX(M)^(CAm-&LOG(XX(m)  )  )/<  (  1  *  0 -B  c  T  A  *  (  G  A  M  -  A  L  P G  (  X 

*  X(M))))*P.P) 

CELT ( M)= ( H/6. C ) *( TP( N,1 )+{^.0*TP(Mt?) ) + ( Z, 0*TP(N, 3 ) ) +T 

*  P  (  N  ,  4 )  ) 

no  2  0-2, 1000 
R(M,i  )=R(N-1,1  )+H 
R ( N, ?)sR( N, I ) +H/2.0 
R(N,3)— R(N,1  )  +H/p,  0 
R  (  N  ,4  )  =R  (  N,  1  )  +H 
T  (  N)  =T(  N-!  )  4-0EL.T  (  N-!  ) 

TPST=T(N ) 

!P(TEST„GC.TZ)  GO  th  9? 

no  2  M= l , 4 

XX ( M)=! , 0~RP*R (M,M) 

7  T  o  (  N  ,  M )  »  X  X  (  M  )  (  G  A M-  AL  0  G  (  X  X  (  M )  )  )  /  (  (  1  *  0  -B  p  T  A  *  ( G  A  M -  A  L  n  G  (  X 

*  XVMM  n*RP) 

2  0Et.T(M)  =  (H/6*0)*CTPfN»l)  +  (  2.0*TP(N,2)  )  +  (  2.0*TP(N,3)  )+T 

*  p ( N ,4 )  ) 

GO  TP  71 

7  p  p  R  =P  (  A’—  1,1) 

TPR  =  T ( M-1  ) 

H= 0.00000! 

N=  1 

P ( N, ! )=R0B 
R  (  N  ,  7  )  ■=  R  (  N ,  !  )  4-H  /  2  •  0 
p  (  m  ,  7  \  -  Q,  (  N ,  T  )  +H/ ?  •  0 
p ( n , 4 ) =p ( N , 1 ) f H 
T ( M ) =TPR 
00  c;  M=!  ,  4 
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t  '  '  r  -  V  - 
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for  nor  ,n*M 

■  f  ? ,  yn  . 

■  '  .  i  '  *  ' 

,  ••  V  *  -  ,,»/)'»  »  .  ‘  ' 

I  ,  )  '  ,  ' 

-  f  ' 

A  ,  f  =  * 


BO 


XX  (  M  )-  1  «,  0-RP*B  (Mf  m  ) 

*>  Tp  (  w,M)  =XX(  M  )  *  (GAm-ALOG(  XX  C  M)|  )  /(  H  .  0-RETA*{G AM-ALO0 (  X 

*  X(M))) }*pD) 

ppLT (^^-(h/6.0)*(TP(N» 1 )+9*0*TP(N,2)+?*0*TP(N,3)+TP(N, 

*  4  )  ) 

PH  IS  N=2,1000 
R  (M,  M=P  (  N-l  , 1  )  +H 
P<N,??=R<Nf 1 ) +H/2.0 
o (N,7  V=R ( M, i  ) +H/P.0 
P(V,4)=R(M,1 ) +  H 

T  t  N) = t ( m - T )+n^LT(N-1  ) 

TFST=T( N ) 

TFf rpcT.Gr.T7)  GC  Tn  70 

nn  t  Ms] t  a 

XX f M  1  =  ]  , 0-R  P*R (M,V) 

7  TP  CN,M) =XX(M) *(GAm-ALOG( XX ( ) ) /< < 1 .6-RETA*<GAM-AL0G( X 

*  x(Min  )*rp) 

6  DBLT  CM)  =  (H/6.0)«( tpc^,i )+2.0*TP(N,2) +2,0*TP(Nf 3)+TP(N, 

*  4  ]  ) 

GO  tp  n 

70  DEI  7  =  ( 1 ,0-pp*o (M, 1  )  ) 

RETURN 

7]  nR.7  =  RR 

RFTIJRN 

7  7  OR  7-1.0 

R F TUP  N 
END 
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FUNCTTTN  P|  N  P  ( R  L  A  M  P ,  K  ) 

DT  =  \1  U'  SQPA 

K  K  =  K  - 1 

P  K=  P  t.nAT  f  KK  ) 

PtMTN=STN(  (RK-M  .0  +  0*  5  >*PI) 

X  =  ( ?.0**0#  666667) /  (3, 0**0*  8333??  ) 

PLN°=  (  °l .  M!N*X*(  PL  AMO**0*3?P  )  /GA  MM  A  (  1 ,  ?333?3  ) 
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SUBROUTINE  VAUJc(7CPIT) 

r  n MM p  T nT  f  TF  ,  TC  f  It  t  C  T  r  F  t  r.  H*>n ,  B  T 

rn^MPM  RL  AM  (40)  ,P.LAMW  (  40)  fCH(  40)  f  CC-C40*  »TK(  40)  f  P.LR(  100 
*  1  ) 

COMMON  R  TP ( 10! ,40 ) ,  PRF$( 101 ) ,  VA( 101 ) 

IF  (  Tl_  .  FQ.  TF  )  GO  TO  0 
TF(7CPTT.lT.1.0F-B)  GO  th  9 
tHFTAF=TC-tF 
THCT1T=TL-Tf, 
nn  3  K  =  i f } 00 
SUM=0.0 
r>n  1  J  - 1 , 40 
P  L  A W - R  l  AMW(  J  ) 
cpp-QL  AW*PL AW*7rR  TT 
IF (CRP*GF.1 74,673 )  GO  Tn  *0 
GO  to  *51 
GO  CHI °=0 • 0 
GO  TO  5  0 

F 1  CHTP  =  r.K(  J  )*RTB(Kf  J)*FXP(-FPP) 

ACHP=ABS (CHI P) 

! F ( ACHP« LT • 1 rOF— ?0)  ^0  5* 

SR  su^=SUM+CHIP 
1  CONTINUE 
5*  V  A  L  U  F 1 =  SU  V 

VA(  K  )=  (  TM^TAF  +  thpt  A!*VALUF  1  )  /  (  T  H  c  T  A  F  T  H  P  T  A  T  ) 

3  r  pnjt  tn'UF 

V  A ( 1 Qi  )=0,0 

pc turn 

3  CONTINUE 

nn  4  K  = 1 , 1 0 1 

V  A  (  K  )  =  1  •  0 
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SUBPPMT INF  Po-SSi ( RT 07 tT7 ) 
r  tmmpm  RT  pt,  TFf  TCf  T !  ,  C  I  C  F  ,  r  H^O »  P  T 

rnMMpM  pL  AM(40)  ,  R  t  A  MW  (  40  )  *  CH  (  40  )  ,0r  (  40)  ,0K(40)  ,Rl  R{  *>00 
*  1  ) 

p HMvnw  RTR  <j  Cl  ,40),RRFSn01),VA(101) 

0!  wcjgc  T  ON  TP  (  !  01  ,  A)  f  n<=LT(  i  n  i  )  t  7 (1 01 , A ) 

H=B  107/100,0 

N=1 

7 ( n, 1 >=0. n 

£<Nf  2)*Z(Nt  1 )  +H/7.0 

Z(N,R)=7(N,l)+H/?.0 

Z(Mf4)=Z(N,l)+H 

PPPC (M) =0,0 
HP  1  Ms] , 4 
X  =  7  (N,M) 

r>Fl_PR  =  PFL7  (X  »T7  ) 

Hp!_  S00  =  nFLPR*0F!  PR 

TP  (  N  t  m  )  -  1  •  0  /  (  D F L  $ Q0 l  SQO  ) 

1  C0NTINUF 

P  ELT  (*’ )  =  ( H/6. 0  )  *  { TP  (  N  *  l )  +7  „  0 *TP (N*2)-*-?,Q*TP(N,R)+tp(N, 

*  4M 

on  ?  m=?,iot 

7.  (  N  » 1  )=Z<N-lfl)+H 

Z(N*?)*7(Nt i ) +H/2.0 

Z(N* 3>=7( N,1 ) +H/P • 0 

7(Mf4)=7(KM )+H 

PRFS(N)  =  PR!r$( N— 1  )  +0FLKN-?  ) 

no  R  M  =  I  f  4 

X  =  7 c  Nf  M) 

0FLPR=0FL7<  X,T? ) 

HPL  SOn  =  nFtPR  *0FLPR 

TP ( N , m ) =1 .0 / ( OFL S Q° SQn > 

R  0  PfSjT  TMI  JR 

PELT  (N)  =  (H/6.0)*( TP(N,' ) +2 . 0*TP ( N, 2 ) ♦? . 0*Tp < N, 3 ) +TP ( N , 

*  4)  ) 

2  CONTINUE 

QCTIJRM 

FNO 


1  f 

•  ,  v  <•  1  |  ’ 

T .  yr  . 

r*i  r  r, 

r  ,  < 

y 

t  . 

\  '  r 

f  ,  f  n  •  «  ,  f  1  ,  '  ,  '  r  '  1 


I  '  ,  M  T 

.  ’  \  .  f  .  M  ’  t  r  ,  ,  ' 

-  ;  n  , .  ,  •  f  *  )  <  t  f  •  M  1 
>  i  ♦  *  '  ,  M"lif 
/  ;(  ‘  )  '' 

'  r  -  »■ 

r* ,  ■  I'  fV 

(  '  T  ,  V  >  i  |  !  '  ~ 

f  n  I  v  n  p  -  pn 

I  ?  ■  '  '  -  )  )  r  p  /  ,  "  -  f  *  .  >  r'  T 

»  j  *  ‘  1  ’  •  1  f  ‘  \  ^  * 

I  (  ' 

r  T  ,  -  1  r 

V  ?  -  .  r  -  M  T  -  (  r  f  IM  r 

r  >  ♦  f  ,  •  '  *  :  I  -  .  »  ' 

I  *  ,  M  t  I  ,  M 

-v 

(  f  _  i  iv  .  r  -  .  • 

'  «  r  ‘  ‘  n 

i  ’  ,  •  >  t  -  y 

(  V  ▼  f  V  |  p  <”  >  -  ’ 

n  i  •  >  s-  •  f  n 

)  \  f  f  ,  '  ' 

*  I  f  ',  f  -  ,  M  '  .  +1  r,  ''  T)  f  ,  \P'-:!  )  '  P 

r  f  a 

'"I  '  T  t  ,40^ 


SURSOHTIMP  p  i  m  (  t?(_  i\Mn  ) 

ni^c^cTHN  p  (  1  00?  ,  4)  ,Pon00?f41fRQP(100?,4)fPPPP(100?,4 

*  )  *  OFt  R  (  TOO?  ) 

i  y  n c q  p  (  i  qo?  ) 

CnvMPN  RIOT,  rpf  TC»  TL  f  C  TCP,CH20f  BT 

rnMMPM  Pt  AM (40) f PL4Mwl40)tCH(4O) *CC<40)  ,CK(40) »RLR( 100 

*  1  1 

rnMMPM  ptr ( i c1 ,4o i , prfs ( ioi 1 1 va n oi i 

H  =  0 . 00 1 0  0  0 
M  =  1. 

R ( N , 1  ) =0*  0 

0  R { M  ,  1 1  =  1.0 

RRP  (KM  1=0.0 

RRPP(M,  i  1=0,0 

R(N,p)=R(N,i 1+H/R.O 

DP  (  M  ,  p  )  =  r p  (  M 1 1  1  +  P  p p  ( Kj ,  1  1  *H /  2  •  0 

PRPCM,?)=PPP(N,i ) +RRPP (N,l ) *H/? . 0 

PR  P  o  (  nj  f  2  )  =  (  P  L  A  MO*  RL  &M0*P  ( N  ,  ?  1  *  (  R  (  N ,  7  )  *R  (  M ,  2  )  - 1  •  0 )  *R  P  {  N 

*  ,  ? ) -PR  P (N, 7  1 

1 )  /  R  (.  N  *  ?  1 

P1N,?1~P(N,1 1+H/P.O 

Pn ( N , ° 1 =RP ( N ,1 1 +PRP (N, 2) *M/?.Q 

PRPfN,?1=°RP(N,l )+RPPP(N,?) *H  /  ?  • 0 

PRPP  (N  ,  •>  )-  (PL  AMD*Rl  AMD*R  (N«,  ^  )  *(  P  (N,  ?  )  *R  (  N,  3  )-1 . 0)  *RP  (  N 

*  , ^1-PPp(Nt?1 
M /d (M,^1 

R(N,41=P(N,1 ) +H 

qp  ( M  ,  4  )  =  RP  (  N  ,  1.  )  +  P R P  (  N  ,  ?  1  *H 

P  R p  ( IS] ,  4  )  =  p  R  P  { M ,  1  )  +PP  p  P  ( M  ,  ^  1  *H 

pd  PP  (  M  ,  4  1  =  (  r  I  AMD*PL  AM[)*p  (N  ,  41  *(  R  (  N,  4}  *R  (  N,  4  )  -1  .  o  1  *PR  (  N 

*  ,  4  1  — R  R  P  {  M  ,  4  ) 

1 1/R ( N , 4 1 

net  r  (M  1  =  (H/A.01  *(RPD  (  Mf  i  )  +?„n*PPp{Nl,?  l+^.O+RRPtN,  *  )  +RD 

*  P  (  N  ,  4 1  ) 

nFLRP(N)=(H/6.01*(PRPP(M,i)+?.0*RRPP(N,?)+?.0*RRPP(N,^ 

*  )+Q?pD(Mf 41 ) 

o  L  P  (  N  1  =  R  R  (  N  ,  1  1 

nn  i  m  =  ->  1 100 1 

R  (M,  1  1  =p  (  N-1  ,1  1+H 

PR ( M  » 1  1  =  RP ( N-T , 1  l+D^LP  CN-1  1 

PPD(  N,  1.)  =t?t? p  (  N-l  ,  1 1  +  DELR  P(  N-1  1 

pppp(N',1  1  =  (  P  L  AMH*RL  AMO*P  ( M  ,  1  1  *  (  P  (N  f  1  )  *R  (  Nf  1  1  -1  •  0  )  *RP  (  N 

*  , M -PRP<N, 1 1 
n/p(N,i ) 

R ( N , p 1 =  P ( N ♦ 1  1 +H/P.0 

RP.  (M*?  1  =  RR  1  N#1  )  +  RRPCN#  1)  *H/?,0 

rrP(M,?)=PRP(N,1 )+PRPP(N,l )  *  H  /  7  •  0 

ppop(Myp)  =  (Rt  \ MO*  RL  ft M0*R1N, ?1*(R<M,?)*R (Nf? 1-1.0 )*RR(N 

*  f -> 1-Rpp (M, ?) 

1  1  /  P  (  N ,  ?  1 

o  <Mf o )  =  P ( m, 1  1 +H/P.0 

pp (Mf^)=RP(Ntl H-RRP<Nf2)*H/?.0 
ppp  (M  ,  ** )  =RRP  (N,  1  )  +  P P  P  P  ( N  ,  *?  )  *H  /  7 . 0 


f  '  •  '  1  \  I  l'»l-  I 

. 

(  ?  '  *  f  1 
I  ~r  r  r  \  >  i  ,  ' 
T  ,  ,  ■  ,  i T  f  ^  !  ,  T  .  T 

f )  i  ,(0- ni‘\!  n  *  f  >  -  '  1  ^  f  c  ' 1  '  i  ,  n  ')  ' A  «*  ’*  •* r 


1  r 

(  T  o  m  v  ,  f  ?  n  ? 

.  '  r  * 1  ) 
o  ,  f  -  (  r  ,  1 

,  f  ■  ,  1  r 

.v  ,^,'n  ' '  .  ■  ' 

r  ,  «  \  f  .  '  *  (  r  f  M  t  ,  M 

.  -  \  '  (  ’  ,  >  r-  •  4  I  r  ,  1 '  '  *  ♦  '  '  " 

|  )  1-1  !  ■  -  f  ,  «  '  'rr  ,  »  f  ’  ,  )  ■  I  f  <  '  '  ,  1  ’ 

'  1  I  , 

(  #  '  1  \  (  r 

.  WU  '  '  ,  *  ’  -  f  ‘  ,  I 

,  \  •  '  •*  .  M  *  '  '  ,  )  '  f  »  ' 

f  ,  r  \  >  (  ,  )  -»  f  r  ,  •  P  :  q  -  f  "  .  4  '  Hi' 

'  )  o  '  f  ■  ,  r  -  (  r  ,  * )  n-  (  r  t  - )  ■>  >  ■  t  t  .  )  o  -  nr  *  |  r-  •  M  A  t  )r|f  ,i»)nn<io 

f  ~ ,  '  c;  -  '  . 

*  '  r.  *  '  f  •'  .  ' 

(  ,  \  Mf,"'  '  ,  ' 

i  '  .  >  ♦  f  ,  M  f  *V  ,  M 

I  \  (r  ,  r.|  ,  !  '  i  ,  )  ,  !  *  i  *  I  )  t  ,  •  \  f  ‘ 

t  >.  f  \  f 

/  .  )  \  '  ■ 

4  I  )  4-  I  ■  '  (  r  )  M  (  '  >  t  ' 

► 

I  '  *  1  * 

,  (  ,  ,  M  t  )  ‘  )  (  ,  v  '  '  ’  1 

'  (  .  "  -  n  *  f 

I  r  f  »  -  •  >  I 

r  r  -  <  r 

t  ,  f  -  • »  '  ;  f  1 : 

f  r -  i  r1  >  f  '  *  f  -  -  (  f ,  '  ' 

f  '  -  )  <  j  (  *  (  r  ,  r  -  '  p  f  '  ,  •  ' 

•  (  I  #  '  M  ,  ’  ^  .  I  '  <  *  , M )  ’  *  1  ‘  *  ’  '  ■  *  r  .  * ' 


i  •  ,  M  \  l 

.  '  \  i  ♦  I  r  •  *  "  *  » 

\  '  f  ,  M  t-  (  r  .  "  >  f  .  )<)1 

r  '  (  r  f  * )  -»  (  *  ,  '  '  -  f r  f  )  ■  °  i 

.  ■  ft  I  '  .  ’ 

('  t  '4  >  '  •  i  . 

f  <•  .  >  i .  v  f  f 


\ '  *  i  ’  ,  1 

r  ,  \  (  t  ,  )  f  r  .  )  -.I’,'  1 


RRP*5  <N,  3  )  =  <  R  t.AMn*Ri.  AMD*R  <N,3)*(R(Nt3)*R(N,3)-l«0)*RR(N 

*  t3)-RRP(N,->) 

1 )/P (N,3) 

R(N*4)=R(Ntl ) +H 

PR  CN,4  )=RR(  N,  1.  )  +RPP  (  M  ,  *  )  *H 
PRP(Mt4)=RRP (M* 1 )+RRPP(N,^)*H 

RRP P  ( N  ,  4  )  =  {  R  l  AMn*P}_  4 MH*R  ( N ,  4 )  *  {  R  (  N  ,  A )  *R  (  N,  4  )  - 1 . 0 )  *° R  (  N 
«  f4)-pPP(M,4) 

1)  /R(N’,A) 

OELR  (N  )  =  (  H/6. 0)  *(  PPP  (N,  1 )  +  2  •  0*PRP  (  N  ♦  ?  l+^.O^PPPIN^l+PR 

*  P( N,4)  ) 

DFLPP(M )-(H/6 .0 ) *( RRPP(N»1 ) +2.0*RPPP(N, 2 ) +  2 . 0*RRP? ( N , 3 

*  |+DOOP(M,4)) 

RLR  {  N)  =°R(N, 1  ) 

1  CONTINUF 
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SU*T?^UT  inf  R  PF  {  7 , G  AM , rft A , RS , RP  ) 
rnviMPM  R  r  nr  ,  tf  ,  Tf  ,  it  t  C  T  C  p  ,  CH^O,  R  I 

C  f)  mv  pm  RL  AM  (40)  ,RLAMW(40)  ,CHf  40)  ,cr{  40)  ,CK(  40)  ,R!.R  non 

*  i  ) 

COMMON  RTP ( 101 ,40) ,PRES (101 ) ,VA (101 ) 

Ic ( 7.PQ.0.0 )  GO  TO  10 
CHT=0.O 
nn  7 4  Kssl  ,40 
PL AM n= PL  AM ( K  ) 

FP  P  =  R  L  A  M 0 * P  I.  AMD*  7 
IF(FRP.GF.174,673)  GO  TO  ^0 
GO  TO  R7 
SO  CH  AP  S  =  0. 0 
GO  TP  c;? 

^7  CHAPSsCH(K)*FXP(-Fpp) 

AC PTK=ABS( CHAPS) 

TF( ACPTK.LT. 1.0F-20)  GO  Tn  ss 
F 2  C H T  =  CHI  + CHAPS 
74  C0NT T NIJF 

ff  TP(tL.FQ.tf)  GO  TO  qi 
TT  =  ( Tp-Tf  )  / (TL-TF) 

GAM=1.0/B T 

rr=(TT*CICF)/(2.0*CH20*CHI  ) 

BET  A  =  1 .0/PR 
X=G AV-BR 
T  F  r  X  )  11,12,12 
11  V  =  A  R  S  (  X  ) 

I F  ( Y  •  G  F  .  1  74  •  673  )  GO  TO  F R 
RS  —  EXPfX  ) 

111  CONTINUE 
PP=1 , 0— R  S 
GO  to  10 
8!  R  S  =  0 . 0 

GAM=1. 0/RT 
RFTA=0.0 
GO  TO  1 1 1 

F  P  R  S  =  0 . 0 

GO  TC  11.1. 

1?  PS»1.0 

WR1tf(6,100)  X 
GO  TO  111 
13  p  S = 1 • o 
GO  Tr  ill 
1 o  RFTyRN 
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SUBROUTINE  FUNC ( X  t GOF  X  , GPp I mX ) 

GOM^GN  BTOTfTFtTC»TLtCICF*CH2ntRT 
PT  =  \U1  B9?.6 

CDNST=  (  (  3,0**0*  666667 ) *GAMMiM  1. 333RR3 )  )  /  (  (  ?.  0**0.  33  3^3 
*  ?3)*GAMMA(0. 

R666667  )  ) 

$T-STN( <X*PI /4.0I-I PI/3.0) ) 

S?=STN( (X*PT /4.0) -(?.0*PT/3.0) ) 

ci=cnsnx*Pi  /4.0)-(r i/3#o n 

C  ? = C  0  S ( ( X*PT /4,0)-( ?.0*PI/7.0) ) 

GOFX=(1.0/BIOT)+( (CONST^S? ) / ( S! * ( X**0. 666666* ) ) ) 

R TMX=( f  S  l*r ?)- ( S2*G1 ) ) /( SI *S1 ) 

Y=-(  ?.0/3.0)*<S?/(Si  *(X**i  .  666667)  )  > 

7  =  ( PI/4*0)*P tmx/ ( X**0.  6666667) 

GPRT  MX  =  COM<;T*(  Y  +  7  ) 
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SUBROUTINE  T R  ANSI  (X,  EPSLON,  M  ) 
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r  . ROOT  FINDER  FOR  TR  ANSCFNDENTAL  FOUATIONS 
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C  . US  AGP 
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f  CAl  t.  TRANS  1  (ROOT  t  ACCRCY  f  TTFRNS ) 


FUNC  TS  A  1JSFR  SUPPt  I  pO  SUBROUTINE  THAT 
CALCULATES  F(X )  AMO  F»(X). 

ENTRY  TR ANSI 


FIND  A  REAL  ROOT  nf  tHp  sc  AL-VALU^D 
TRANSCENDENTAL  FQUAT  T  ON  f(X)=0. 
ROr'T»,,,,.!NTTTAL  GUFSS  AND  CALCULATED 
C  7  ERO 

C  ACCRCY. . ..ITERATION  BOUND 

C  XU  +  n  -X  (T  XEPSLON 

f  TTFRNS, .. .MAXI  MUM  NUMBER  nF  ITcrattoms 
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XPP=X 

CALL  FUNC( XPP»YPP, YDPP) 

Xp  =  X pn-YPp/YOPP 

nn  o  f= 1 , M 

PALL  FUNCf XP,YP,YQP> 

D=XP-XPP 

!F( A8S(D) .LT.2.98E-08)  GO  Tn  ?4 
A=  <  P. 0*Y0P  +  YDPP-a.0*( YP-YPP ) / D )  / D 
IJ  =  Y  P  /  Y  0  P 

y=xD-u* ( 1 .0+U*A/YDP) 

IF  (  ABS  MX'XP)  /  X  J-EPSLON)  T,  B  ,4 

4  x  PP  =  X  p 

XP  =  X 

YopryP 
2  yopp=yhp 
WRTTF(6,'’'7) 

Rp  format  {  *  fAILURF  —  INCREASE  w  AND/OR  INCREASE  FPSIONM 
f  q  t  n  P 

24  W  R I T F ( 6  »  ° 3 ) 

23  format ( «  FAILURE — ACCURACY  REQUIREMENTS  EXCEED  S • p •  S I 
*  GN  T  f IC  ANC F  * ) 
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SUBROUTINE  MJWFR4 ( A,B ♦  SS4 ) 

COMMON  RIOT  tTFtTC»TLfC!CF*CH20»BT 

CHMMfN  PL  AM ( 40 ) ,RLAmw(40) ,CH{40) ,CC(40) ,CK(40) ,RLR< TOO 

*  l ) 

COMMON  RTR.1 101  f  40)  ,PRES(101  )  tVM  101 ) 

F( r )=r*( i „0-C*C ) 

on 

N  3  =  N  n  /  ? 

H=(  R-A  ) /FLOAT  {  NO  ) 

M1=IFIX< A*1000.0)+l 

mp=tptX( { A+H) *1000.0 )+l 

S=F  (  A)  *RLR  <  Ml  )  +4,0*F  (  A  +H  )  *R  LR  {  M  ?  ) 

NT =NR-l 
OH  1  T=TtNl 

m  0=  T  F  T  X  (  {  A+H  * FLn  A T { 9  *  t  >  )*1  000*0+1 
M4=icyx(  ( A+H* FLOAT!?* I +1))  *10000)+! 

c  =  9, Q*r( A+H* FLOAT ( ?*I  )  ) *R|.  R ( mR ) +4* 0*F ( A+H* FLO  AT (?*I+1) 

*  )*RLRCM4)+$ 
i  CONTTNUF 

VR=TFIX(B*1 000.0)+! 

X=f(B)*RLR(m*S)+F 

SS4=H*X/*.0 

RETURN 

END 
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roMMOisj  RTR  ( 101 ,40)  ,PRES(  101  )  ,VA(101 ) 

F  (  0  )  =  0  *  { 1 . 0-0*0 ) 

mr=i on 
NO  =  M?/'> 

H=( s-A ) /FLPAT(N? ) 

Ml  =  T  F I X ( A*)  000.0) +1 
M?=!PTX( ( A+H ) *1000.0 ) + l 

S= F ( A  )  *  R  L  R (Ml )*RLR (Ml ) +4.0*F( A  +  H)*Rl_R (  M2)*RLR(M2) 
m 1 =M?- 1 
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mo=TFTX(  (  *  +  H*FLPAT(?*T  U*1  000.0)+! 

M4=  T p I  X ( ( A+H* FLOAT (?*!  +  l) ) *1000.0)+! 

<;=?,  0*F  (  A+H*P(,0  AT  (  )  )  *rlr  (M3  )  *RLR  (  MB)  +4.0*F  (  A+H*FL0A 

*  T( 2*1  +  1 ) ) *R  L 
1R (M4) *PLR(M4)+R 
1  FONT  T MljP 

MF=TFTX( B  * 1 000.0) +1 

X  =  F  (  B  )  *R(_R  (  M5  )*J?LR  (  MB  )  +S 
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SUBPHUTTNF  TNTPG2(A,RtSS?) 
r  pmvpm  Rinr,  TPf  tc,  TL,  r.TCF^CH^n,  BI 

fPMMOM  R|  AM(40)  ,RLAMW(40)  tCH(40)  ,cc.(40)  ,CK(40)  tPLR(  100 
*  1  ) 

COMMON]  RTR  (Id  ,40) ,PRFS(101 ) ,VA(101  ) 

F ( C ) ~C  * ( 1  •  0-C  *C ) 
m  o  =  l  0  0 

N  3  =  N  ?  /  R 

H=( Q- A )/*=L0AT(Np ) 

Ml = TFT X( A*1 000. 0 ) +1 
N1=TF  !X( A*! 00.0)+) 

M'’=  T  F  T  x  (  f  A  +  H  )  *  1  000*0)  +  ! 

N2=IF  !X(  (A+H  ) *1  00,0)+! 

R=F ( A ) *RLP ( Ml )*VA(N! ) +4. 0**= ( A+H ) *Rl_P ( M? ) *VA ( N? ) 

N 1  =  M3- 1 

np  i  t=i,N! 

M  3  =  T  R  T  X  (  ( A+H*FLOAT( ?*I ) ) *1  000*0 ) +1 
N3=  TFI  X(  (A+H*FLOAT  (2*1  H*100«0)+I 
M4=  T  F  T  X  (  (  A  +  H«FLOAT( '>*!+-!  )  )  *1  000*0)  +  ! 

N!4= T  F  T  X (  ( A  +  H*FLOAT( ?*I  +  ! ) ) *  1 00. 0 )  +  1 

F  =  3. 0*F(  A+H^FLHAT ( ?*T  ) ) *PLD (M3 ) *VA ( M3) +4. 0*F { A+H*Fl 0*T 
*  ( 2*1+1 )) *RLR 

1  (  M4  )  *\/ A  (  M4  )  +  $ 

1  CDNTTNUF 

mF=Tctx( 0*1 000*0)+! 

=  T  F  T  X ( R*1 00*0)  +  ! 

X=P(B)*RL°(M6)*VA(N3)+S 

RS'>-H*X/3.0 
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APPENDIX  B 

ORDER  OF  MAGNITUDE  ANALYSIS  OF  THE  GOVERNING  EQUATIONS 

Consider  the  governing  differential  equations:  continuity, 
momentum,  energy.  For  cylindrical  coordinates  assuming  axi -symmetry , 
no  sources,  and  negligible  viscous  energy  dissipation. 

Continuity: 


h£+p&  +  f!fK>]  =  0 


B-l 


Momentum: 


3v  3v  3v 

_ _  +  w  — —  +  v  — — 

3t  z  3z  vr  3r 


1  Sn  3^,  1  3V7 

Jlf+v[r^  +  F  — +  — ]  -  9 


3r 


3z 
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+  V 


+  V. 


=  .1M 


32v 


z  3z  r  3r 


P  3r  +  vCZ2 
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n  3V 

r  1  r 
—  +  - - -  + 


r  3r 
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Energy: 


329  ,  1  39  ,  3^i  =  1  rli  +  »  99  +  ■■  99 

2  r  3r  „  2  a  L3t  r  3r 
3r  9z 


vz  3l] 
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Consider  the  energy  equation  B-4.  Normalizing  we  obtain 


.6  v  +ii  v  +  ££  1  .  fl  +  4i  +  fi 
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[\][^f  +  1  |f]  +  [\]  =  [!|_]  |i 

R/  9R2  K  9K  Z  2  9Z2  atc  9t 

c  c 

e  u  e  v 

+  r— c-  C1  u  r  c  ci  v  ■£& 

LaR^  J  U  9R  LaZ  J  V  9Z 


where 


e 


V 


V 


^  0  5  ^  D  »  Z  -  7  »T  4-  >  U  -  1 1  »  V  =  w 


R 


U 


V 


Multiplying  the  above  equation  by  — U  yields 

R  2 
c 


a2.  -i  rv ,  R  2  ,2 .  R/  _  .  UR  ~ ,  V  R 
9  x  +  1  iA  +  r  JLi  LA  =  r_9-  .i  iA  +  r„c.  c.~|  n  Lk  +  r  c  c  ~i  v  M 

3R2  R  9R  +  LZcJ  9Z2  LatcJ  9t  L  a  J  U  9R  L  7  J  V  "7 


9Z 


or 


R  2  ,2 


R 


9R2  R  9R  LZc  9Z 


•at  9t 
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LA  +  1 LL  +  r— ' l*"  LA  =  r'.A..i  JA  +  r.AR.c.  ~i{y  A  +  rlLAi  n  A) 

9  T  d  L7  -I  2  Lat  J  9t  L  aZ  Jlv  97  L\/  d  J  u 


U  Z 


9Z  LV  R 
c  c 


9R- 


at 

If  — j  »1  ,  then  the  unsteady  term  is  insignificant.  Therefore 
Rc 


92({)  1  9<J>  _  r^C  c  -1  ry  9(j)  ,  r  C  C-I  M  9(j)  x 
^2  +  R  9R  -  L  aZc  J  {V  9Z  +  LVcRcJ  U  9R  1 


where  R-  «  Z  ,  i.e.  the  axial  conduction  is  negligible.  If  advection 
c  c 

2 

is  as  important  as  the  conduction  terms,  then  VcRc  /aZQ  =  0(1).  Hence 


2 

R  v 

z  =  S — £  =  r  pr  ReD 
c  a  c  R, 
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The  energy  equation  becomes 


iLiL  +  J_M-=  {y  f  c  C1  n  li}  B_5 

.d2  R  9R  1V  9Z  LV  R  J  U  9RK  b  b 

d  K  C  C 

The  solution  of  equation  B-5  is  valid  for  Z  of  0(1);  however,  the 
solution  is  expected  to  be  accurate  for  smaller  values  of  Z  and  therefore 
any  solution  will  be  extrapolated  into  this  region. 

Now  consider  the  continuity  equation  but  rewritten  in  the  form 


9p 

9t 


rrr  +  V. 


+ 

9r 


9£ 

9z 


p[ 


9v. 

dz 


if-  (rv  )]  = 
r  9r  r ' 


0 


Normalizing  we  obtain 


where 


But  if 


then 


-  ^cPcP' 


9R 


and 


-  3ecpcp' 


9Z 


where  $  is  the  coefficient  of  thermal  expansion 


•  o  *  [(/•''  +\^  sv  ’  :j  V 


.  Mh  .q  [2  :  ‘q  [2-f]  > 
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and  the  normalized  continuity  equation  becomes 


r-k-1  Ml  +  D.  +  r^i  ^RUh 

LtcVcJ  3x  +  p  13Z  LVcRcJ  9R  J 

-  69c {v  If +  u  i}  • 

c  c 


For  the  unsteady  term  to  be  insignificant  , 


B-6 


»  1 


as  in  the  energy  equation.  For  water  under  atmospheric  pressure 
(36  )  _  «  1  and  the  right  hand  side  of  equation  B-6  is  suppressed 

C  tTlaX 

since  the  terms  inside  the  brackets  {}  are  of  order  0(1)  as  shown  by 
the  analysis  of  the  energy  equation.  Since  p'  =  0(1) 


3V  .  rUcZCi  3 ( RU )  _  n 

3Z  LV  R  J  3R 
c  c 

u  z  Rc 

Taking  =  0(1)  we  define  Uc  =  (— )  Vc  ,  therefore 

CC  C 


9V  ,  3 (RU)  _  n 
3Z  9R 
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Next  consider  the  momentum  equations 


.  *»mng»»nt  *1  o*  ml  .*•**»“  •« 
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let 


p(r,z)  =  p°(z)  +  p'(r,z) 


where  p°  =  hydrostatic  pressure  at  an  interface 


p  =  pressure  departure  from  p 


Thus 


3z 


9£1+  1L 

3z  3z 


or 


3£ 

3z 


.pg+i 

Mwy  3z 


The  axial  momentum  equation  B-2  becomes 


3v,  3v  i  ~ 

— -  +  v  — -  +  v  — —  =  -  — 

3t  z  3z  r  3r  p  3z 


32v  3v_  32v_  pi( 

3r  K 


B-8 


By  definition,  the  coefficient  of  thermal  expansion  is 


3  = 


1  (i£) 

p  '8T 


or 


g  -  1 

B  p  Tw  ■  T  • 


Therefore  equation  B-8  becomes 


.  n  q'  a  &  1  *3  1  /  •  *  90  -  3  -  !C  f  -00  ^ 
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3v. 

t 

r  9r" 
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where 


0  =  T  -  T 


w 


Normalizing  the  momentum  equations,  we  obtain 


[^]f +[^-]V§+[^]uf 


vVCn  ,32V  .  1  3V 


A,  32V 


^  If  +  ^  +  iW>  +  ^  ^  + 

*C  c 
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pRc  3R 


R 


3R 
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II  vUr  a^ll 

4) +  [-7]  ^4 

FT  Z  c  DT 

c 


where 


P'  = 


:  JL 


I  • 


R 


Taking  the  axial  momentum  equation  and  multiplying  by  ,  we  obtain 
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c  c 
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Since  the  pressure  and  viscous  forces  are  important 
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T  -  T  *  6 


9i9r1w 
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P  'R/ 

[^Z77]  =  0(1) 

C  C 


defining 


D  ,  yZCVC 

c  =  R2 


If  the  unsteady  term  is  to  be  insignificant 


V 

R  2 


»  1 


where 


t  a  tv  n 

c  _  c  ^  ! 

R  2  R  2  Pr 

c  c 


at 


vt 


For  water  Pr  >  1  ,  therefore  if  — %  »  1  then  — %  »  1.  Hence  if  the 

R  /  R/ 

c  c 

energy  equation  is  quasi-steady ,  so  are  the  momentum  and  continuity 


equations . 


The  axial  momentum  equation  then  becomes 


[• 


2 

Rc  Vc1n,  3V  .  ..  3V-,  _  3P  .  32V  .  1  8V  ,  r  1  -I2  3 2V  .  rGrn  , 
][V  3Z  +  U  3R]  -  '  3Z  +  ^2  +  R  3R  +  [P?R?]  ^2  +  W/ 


vZ 


8Z  3R2 


where 


Gr 


R 


6gRc30c 


V 


,  Grashof  number 


Re 


R 


V  R 
c  c 

V 


,  Reynolds  number  . 


Therefore  if  p-^—  «  1,  i.e.  axial  conduction  negligible,  the  equation 


reduces  to 


o>  q 


0 


o'  o'1 


. 


N 
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[p?]  [v  lr +  u  = 


1E1  +  l_i  4.  1 4.  r5r.n  * 

3Z  R2  R  3R  Re  R  * 

c 


B-9 


Normalizing  the  radial  momentum  equation  B-3,  we  obtain 


[J_]  [_J ]2  [u  QL  +  v  — ]  =  -  — 

LPrJ  LPrReJ  Lu  dR  +  v  3Zj  9R 


+  [_1 — ]2  r^L  +  1  20L 

LPrReJ  l8R2  R  SR 


4  2 

—1  +  [— L_i  iLLL 

R2J  LprReJ  SZ2  * 


Again,  since  -p-R--  «  1  then 


SP 

SR 


=  0  . 


B-l  0 


The  governing  equations  will  now  be  applied  to  the  ice-free 


and  freezing  zones. 


A.  Ice-Free  Zone 


Taking  Rc  =  r^  and  Vc 


=  JL 


TTr 


y  equation  B-5  becomes 


0 


3^  +  i!!i=  v  ^ 


SR 


R  R 


SZ 


UcZc 


remembering  the  flow  is  fully  developed  and  y— ^ —  0(1)  from  continuity 


c  c 


Since  the  advective  terms  are  as  important  as  the  conduction  terms 


r0  Pr  Rero 
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tcao 

The  unsteady  term  in  the  energy  equation  was  neglected  although  — —  _>  1 

r0 

experimentally.  The  tacit  neglect  of  the  unsteady  term  in  the  energy 
equation  reveals  that  the  solution  for  the  temperature  distribution  in 
the  ice-free  zone  is  a  steady  state  approximation. 

Since  the  flow  is  fully  developed  hydrodynamical ly  and  if 
natural  convection  is  negligible  i.e. 


Re  »  Gr 
ro  ro 


where 


Gr 


sgr0  e, 


ro  2 

V 


and 


ec=  -  =  <Tl  -  W'O  -  V’z+)]max 


then  the  momentum  equation  would  reduce  to 

92V  1  3V  _  9 Ps 
3R2  R  3R  '  3Z 


B. 


Freezing  Zone 
1 .  Liquid  Phase 
Taking  RQ  = 


and  V  = 


+ 1 

3R2  R  3R 


2 


the  energy  equation  becomes 


9<J) 


£ 


9Z 


+  U 


9R 


B-ll 
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Consider  the  energy  balance  at  the  ice-water  interface 


ST 


-  k 


£ 


£  Sr 


d^>  . 
i  dt 


Normalization  reveals 


where 


£ 

Lk 

Ki 


a] 


S<|> 

SR“ 


£ 


+  [ 


piLiRc 


ki6ci 


,  dA  _  8(N- 
t  J  dx  ~  SR 


A 


a 


e 


c£ 


0  . 
Cl 


the  superheat  ratio 


B-12 


For  a  growing  layer  of  ice, the  coefficient  of  the  middle  term  of  equation 

B-12  cannot  be  very  small  without  suppressing  the  nature  of  the  problem, 

ice  formation,  therefore  consider  this  coefficient  of  order  0(1).  The 

significance  of  the  heat  flux  in  the  water  depends  on  the  coefficient 

k£  - 1 

-r^~  o  ;  for  the  experimental  cases  this  term  is  of  the  order  (10  ).  This 

Ki 

indicates  that  the  exclusion  of  the  unsteady  term  in  the  energy  equation 
has  only  a  small  effect  on  the  solution  for  the  interface  history. 

Considering  the  axial  momentum  equation,  the  inertia  terms  are 
insignificant  when  Pr  »  1  (Pr  =  13.6  for  freezing  water)  and  natural 
convection  is  negligible  if  RerQ  »  GrrQ 


,  where 
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0 


"  Tf^max  [^(0,z*)  (T 


0 


-  Tf)  -  T-] 
f  fJmax 


2.  Solid  Phase 


Consider  the  energy  equation  as  given  by  equation  B-4;  we 
obtain  after  normalization 


d2(p. 


1  Wl  ,RrJ  d‘ %  A  _ 


9R 


1L  +  1  __1  +  r_£i  Hi  _  r  c  u^i 

2  R  9R  LZcJ  az2  "  [t^7]  97"  • 


c  1 


Consider  Z  equal  to  the  frozen  length  of  the  tube,  i.e.  Z  =  £  ,  and 

c 

Rc  =  rQ  then  axial  conduction  in  the  ice  is  negligible  if  rQ/£  «  1  , 
yielding 


A  i!V  rCo!_T  Zli 

3R2  R  3R  Ltca.J  3x 


B-13 


As  demonstrated  previously,  by  normalizing  the  equation  for  the  energy 
balance  at  the  solid-liquid  interface  an  expression  for  t  is  found, 
such  that  equation  B-13  takes  the  form 


32*.  ,  3<f>,  3*, 

- -  +  - - -  =  ste  — - 

9R2  R  9R  3t 


If  Ste  «  1,  the  unsteady  term  is  negligible  and  the  solution  is  "quasi¬ 
steady"  . 

If  the  ice-water  interface  has  reached  a  steady-state  value 
the  equation  is  simply 

32*i  ,  3^ 

3R2  "R3R 


=  0  . 
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